On the construction of high-order force gradient algorithms 
for integration of motion in classical and quantum systems 
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A consequent approach is proposed to construct symplectic force-gradient algorithms of arbitrar- 
ily high orders in the time step for precise integration of motion in classical and quantum mechanics 
simulations. Within this approach the basic algorithms are first derived up to the eighth order by 
direct decompositions of exponential propagators and further collected using an advanced composi- 
tion scheme to obtain the algorithms of higher orders. Contrary to the scheme by Chin and Kidwell 
[Phys. Rev. E 62, 8746 (2000)], where high-order algorithms are introduced by standard iterations 
of a force-gradient integrator of order four, the present method allows to reduce the total number 
of expensive force and its gradient evaluations to a minimum. At the same time, the precision of 
the integration increases significantly, especially with increasing the order of the generated schemes. 
The algorithms are tested in molecular dynamics and celestial mechanics simulations. It is shown, 
in particular, that the efficiency of the new fourth-order-based algorithms is better approximately 
in factors 5 to 1000 for orders 4 to 12, respectively. The results corresponding to sixth- and eighth- 
order-based composition schemes are also presented up to the sixteenth order. For orders 14 and 
16, such highly precise schemes, at considerably smaller computational costs, allow to reduce un- 
physical deviations in the total energy up in 100 000 times with respect to those of the standard 
fourth-order-based iteration approach. 

Pacs numbers: 02.60.Cb; 05.10.-a; 95.10.Ce; 95.75.Pq 



I. INTRODUCTION 

Understanding the dynamic phenomena in classical 
and quantum many-body systems is of importance for 
the most of areas of physics and chemistry. The de- 
velopment of efficient algorithms for solving the equa- 
tions of motion in such systems should therefore impact 
a lot of fields of fundamental research. During the last 
decade a considerable activity has been directed on 
the construction of symplectic time-reversible algorithms 
that employ decompositions of the evolution operators 
into analytically solvable parts. The decomposition algo- 
rithms exactly preserve all Poincare invariants and, thus, 
are ideal for long-time integration in molecular dynam- 
ics jl(J and astrophysical simulations. The reason is 
that for these algorithms the errors in energy conserva- 
tion appear to be bounded even for relatively large values 
of the size of the time step. This is in a sharp contrast to 
traditional Runge-Kutta and predictor-corrector schemes 
where the numerical uncertainties increase lin- 
early with increasing the integration time |]|l4] |l7j . 

The main attention in previous investigations has been 
devoted to derive different-order decomposition algo- 
rithms involving only force evaluations during the time 
propagation. For instance, the widely used velocity- and 
position- Verlet algorithms [p^|Jig| l relate, in the general 
classification, to a three-stages decomposition scheme of 
the second order with one force evaluation per step. The 
fourth-order algorithm by Forest and Ruth || corre- 
sponds to a scheme with three such force recalculations 



and consists of seven single-exponential stages. Sixth- 
order schemes are reproduced beginning from fif- 
teen stages and seven evaluations of force for each body 
in the system per given time step. With further increas- 
ing the order of force decomposition schemes, the num- 
ber of stages and thus the number of the correspond- 
ing non-linear equations (which are necessary to solve 
numerically to obtain the required time coefficients for 
single-exponential propagations) increases drastically. In 
addition, such equations become too cumbersome and 
all these, taking into account the capabilities of modern 
supercomputers, led to the impossibility of representing 
the direct decomposition algorithms of order eighth and 
higher in an explicit form [5|. In order to simplify this 
problem, it was proposed p], 3 ,p|,p|j2?^-|2^| to derive higher- 
order integrators by composing schemes of lower (actu- 
ally second) orders. The resulting second-order-based 
composition algorithms have been explicitly obtained up 



to the tenth order 
Relatively recently 



3,221. 

24-|2q], a deeper analysis of the 
operator factorization process has shown that the class 
of analytically integrable decomposition integrators can 
be extended including additionally a higher-order com- 
mutator into the single-exponential propagations. As a 
consequence, a set of new so-called force-gradient algo- 
rithms of the fourth order has been introduced. A distin- 
guishable feature of these algorithms is the possibility to 
generate solutions using only positive values for time co- 
efficients during each substage of the integration. This is 
contrary to the original decomposition approach, where 
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beyond second order (as has been rigorously proved by 
Suzuki |§) any scheme expressed in terms of only force 
evaluation must produce some negative time coefficients. 
We mention that applying negative time propagations 
is impossible, in principle, in such important fields as 
non-equilibrium statistical mechanics, quantum statis- 
tics, stochastic dynamics, etc., because one cannot sim- 
ulate diffusion or stochastic processes backward in time 
nor sample configurations with negative temperatures. 
In the case of stochastic dynamics simulations it has 
been demonstrated explicitly |27],^] that using fourth- 
order force-gradient algorithms leads to much superior 
propagation over standard Verlet-based schemes of the 
second order in that it allows much larger time steps 
with no loss of precision. A similar pattern was ob- 
served in classical dynamics simulations comparing the 
usual fourth-order algorithm by Forest and Ruth with 
its force-gradient counterparts |£6|. 

Quite recently, Chin and Kidwell j2jj has considered a 
question of how to iterate the force-gradient algorithms 
to higher order. The iteration was based on Creutz's and 
Gocksch's approach |50| according to which an algorithm 
of order K + 2 can be obtained by triplet construction 
of a self-adjoint (i.e. time-reversible) scheme of order K . 
Then starting from a fourth-order integrator, it has been 
shown in actual celestial mechanics simulations that for 
orders 6, 8, 10, and 12, the numerical errors correspond- 
ing to the force-gradient-based schemes are significantly 
smaller than those of the schemes basing on iterations of 
usual non-gradient algorithms. The resulting efficiency 
of the integration has also increased considerably despite 
an increased computational efforts spent on the calcula- 
tions of force gradients. The same has been seen in the 
case of quantum mechanics simulations when solving the 
time-dependent Schrodinger equation pl| . 

It is worth emphasizing, however, that the iteration 
scheme by Chin and Kidwell is far to be optimal for deriv- 
ing high-order integrators belonging to the force-gradient 
class. The reason is that the number of total force and its 
gradient evaluations increases too rapidly with increas- 
ing K. Remembering that such evaluations constitute 
the most time-consuming part of the calculations, this 
may restrict the region of applicability of force-gradient 
algorithms to relative low orders only. Note that high- 
order computations are especially desirable in problems 
of astrophysical interest, because than one can observe 
over a system during very long times. They may also be 
useful in highly precise molecular dynamics and quantum 
mechanics simulations to identify or confirm very subtle 
effects. 

In the present paper we propose a general approach 
to construction of symplectic force-gradient algorithms 
of arbitrary orders. The approach considers the splitting 
and composing of the evolution operators on the basic 
level, taking into account the explicit structure of trun- 
cation terms at each given order in the time step. This 
has allowed us to obtain exclusively precise and econom- 
ical algorithms with using significantly smaller number 



of single-exponential propagations than that appearing 
within standard decomposition and iteration schemes. 
The paper is organized as follows. The equations of mo- 
tion for classical and quantum systems are presented in 
section II. A. The integration of these equations by di- 
rect decompositions and their force-gradient generaliza- 
tion are described in section II. B. Explicit expressions 
for basic force-gradient algorithms of orders 2, 4, 6, and 
8 are also given there. The higher-order integration bas- 
ing on advanced compositions of lower-order schemes is 
considered in section II. C. The composition constants for 
fourth-, sixth-, and eighth-order-based schemes are cal- 
culated and written down in the same section up to the 
overall order 16. Sections III. A and III.B are devoted 
to applications of obtained force-gradient algorithms to 
molecular dynamics and celestial mechanics simulations, 
respectively. A comparative analysis of the new algo- 
rithms with existing integrators is made there as well. 
The final discussion and concluding remarks are high- 
lighted at the end in section IV. 



II. GENERAL THEORY OF CONSTRUCTION OF 
FORCE-GRADIENT ALGORITHMS 

A. Basic equations of motion for classical and 
quantum systems 

Consider first a classical A^-body system described by 
the Hamiltonian 
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1 N 



T + U, 



(1) 



where is the position of particle i moving with velocity 
Vi = dr^/dt and carrying mass rrii, ip(rij) = y(|r, — ry|) 
denotes the interparticle potential of interaction, and T 
and U relate to the total kinetic and potential energies, 
respectively. Then the equations of motion can be pre- 
sented in the following compact form 



dp 

dt 



= [poH} = Lp(t) . 



(2) 



Here p — {r, v} = {r,, Vj} is the full set (i = 1,2,..., N) 
of phase variables, [o ] represents the Poisson bracket and 
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(3) 



is the Liouville operator with f, = — Y^jCjjti) v' ( r ij) r ij / r ij 
being the force acting on particles due to the interactions. 

In the case of quantum systems, the state evolution can 
be described by the time-dependent Schrodinger equation 



if^±=H{v)^={T 



U{v)U, 



(4) 
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where T = ^|Eti'' 2 ^i/ m i an d U are the kinetic 
and potential energy operators, respectively, and ip is 
the wave function. So-called quantum-classical dynam- 
ics models ||| can also be introduced. This leads to 
a coupled system of Newtonian (2) and Schrodinger (4) 
equations. But, in order to simplify notations, we re- 
strict ourselves to the above purely classic and quantum 
considerations. 

If an initial configuration p(0) or ip(0) is provided, the 
unique solution to Eq. (2) or (4) can be formally cast as 



n(t) = e""ft(0) 



,£Ai 



11(0) 



(5) 



where At and I —t/At are the size of the single time step 
and the total number of steps, respectively, 71 denotes ei- 
ther p or ip, whereas C corresponds to L or —iH/fr. As is 
well known, the time evolution of many-particle systems 
cannot be performed exactly in the general case. Thus, 
the problem arises on evaluating the propagator e £At by 
numerical methods. 



B. Integration by direct decompositions 

1. Original decomposition approach 

The main idea of decomposition integration consists 
in factorization of the full exponential operator e £At on 
such subpropagators which allow to be evaluated analyt- 
ically or at least be presented in quadratures. Within the 
original approach, this is achieved by splitting the oper- 
ator C = A + B into its kinetic A and potential B parts, 
where A = v-d/dr or A — —iT/h and B = a-d/dv with 
a = {a.;} = {fj/mj} being the acceleration or B = —UA/h 
for the cases of classical or quantum mechanics, respec- 
tively. Then, taking into account the smallness of At, 
the total propagator can be decomposed |j]-|||||4j using 
the formula 



^{A+B)At+0(At K + 1 ) 



P 

n 

P =i 



e Aa p At e Bb p At 



(6) 



where the coefficients a p and b p are chosen in such a way 
to provide the highest possible value for K > 1 at a given 
integer number P > 1. As a result, integration (5) can 
performed approximately with the help of Eq. (6) by ne- 
glecting truncation terms 0(At K+1 ). The precision will 
increase with increasing the order K and decreasing the 
size At of the time step. 

As can be verified readily, the exponential subpropa- 
gators e Ar and e Br , appearing in the right-hand-side of 
Eq. (6), are analytically integrable for classical systems. 
Indeed, taking into account the independence of v on r 
and a on v yields 



e Ar p ee e rv-a / ar {r, v} = {r + vr, v}, 
e BT p ee e ra,a / Sv {r, v} = {r, v + ar} 



that represent simple shift operators in position and ve- 
locity spaces, respectively, with r being equal to a p At or 
bpAt. For quantum mechanics propagations, the kinetic 
part e Ar = e~ lT ' r / Ti will require carrying out two, one 
direct and one inverse, spatial Fourier transforms |sT| , 
whereas the calculation of e Br = q- itU / k is trivial. 

In view of decompositions (6), one can reproduce inte- 
grators of various orders in the time step. In particular, 
the well-known second-order (K = 2) velocity- Verlet al- 
gorithm piP] 



(8) 



c (A+B)At+0(At 3 ) = e B^ e AAt e B^ 



(J) 



is readily derived from Eq. (6) by putting P — 2 with 
ai = 0, fri = &2 = 1/2, and a% = 1. The fourth-order 
(K = 4) algorithm by Forest and Ruth (2| is obtained 
from Eq. (6) at P = 4 with a\ = 0, 02 = = 8, 
a 3 = (1 - 26»), 61 = 64 = 6»/2 and b 2 = b 3 = (1 - 0)/2, 
where 9 = 1/(2 — \/2~). Schemes of the sixth order 
(K = 6) are derivable starting from P = 8 with nu- 
merical representation of time coefficients 0,^]. 

The original decomposition approach has, however, a 
set of disadvantages. First of all, it is worth pointing out 
that with further increasing the order of integration (6) 
to K = 8 and higher, the number 2P of unknowns a p 
and bp begins to increase too rapidly. This leads to the 
impossibility of representing algorithms of such a type 
for K > 6 in an explicit form |J , because it becomes im- 
possible to solve the same number of the resulting cum- 
bersome non-linear equations (with respect to a p and b p ) 
even using the capabilities of modern supercomputers. 
Another drawback consists in the fact that for K > 2 it 
is impossible || at any P to derive from Eq. (6) a de- 
composition scheme with the help of only positive time 
coefficients. For example, in the case of Forest-Ruth in- 
tegration, three of eight coefficients, namely, 03, 62, and 
63, are negative. As was mentioned in the introduction, 
schemes with negative time coefficients have a restricted 
region of application and are not acceptable for simulat- 
ing non-equilibrium, quantum statistics, stochastic and 
other important processes. Moreover, for schemes ex- 
pressed in terms of force evaluation only, the main term 
0(At K+1 ) of truncation uncertainties appears to be, as 
a rule, too big, resulting in decreasing the efficiency of 
the computations. 



2. Generalized force- gradient decomposition method 



From the afore said, it is quite desirable to introduce 
a more general approach which is free of the above dis- 
advantages. At the same time, this approach, like the 
original scheme, must be explicit, i.e., lead to analytical 
propagations. In addition, it is expected that the already 
known decomposition algorithms should appear from it 
as particular cases. 
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Let us first analyze the structure of third-order trun- 
cation errors 0(At 3 ) of the velocity- Vcrlet algorithm in 
detail. Expanding both the sides of Eq. (8) into Taylor's 
series with respect to At, one finds 



0{At 3 ) = 



±[A, [A,B]] + [A,B]]} At 3 + 0(At 5 ) 

(9) 



where [ , ] denotes the commutator of two operators. 
Taking into account the explicit expressions for operators 
A and B it can be shown that one of the two third-order 
operators in Eq. (9), namely [B, [A,B\\, is relatively sim- 
ple and, that is more important, it allows to be handled 
explicitly, contrary to the operator [A, [A, B}]. In the case 
of classical systems it can be obtained readily that 



d 



d 



(10) 



where gj Q — 

2 J2jp fjp/rrijdfia/drjp. In view of the ex- 
pression f iQ = - J2j(j^i) <P'( r ij){ria - rja)/rij for forces, 
the required force-gradient evaluations diia/drjp are ex- 
plicitly representable, i.e., 



N 



Si = -2 ^ 



x ' '13 'ij x 

x (r ir (ai - a,)) 



(11) 



= s( r u) = gi( r )' 



As can be seen easily from Eqs. (10) and (11), the op- 
erator C commutes with B = a-d/dv, and, in addition, 
the function G like a does not depend on velocity. Then 
the force-gradient part CAt 3 /24 of truncation uncertain- 
ties (9) can be extracted by transferring them from the 
left-hand-side of Eq. (8) to its right side and further sym- 
metrically collecting with operator B under exponentials. 
This yields the following force-gradient version 



^A+B)At+0(At 3 ) 



48 £ 



e 



(12) 



of the velocity- Ver let integrator, where already 0{At 3 ) = 
[A, [A,B]]At 3 /\2. 

In the case of higher-order (K > 2) integration (6), 
the operator C will enter into truncation uncertainties 
0(At K+1 ) by various combinations. They can be ex- 
tracted similarly as for K = 2, and we come to a force- 
gradient decomposition approach. The most general rep- 
resentation of this approach is 

p 

c (A+B)At+0(At K + 1 ) _ J~J c Aa p At c Bb p At+Cc„At 3 ^g-j 
p=l 

where again at a given P the coefficients well as 

c p must be chosen in such a way to cancel the truncation 



terms 0(At K+1 ) to the highest possible order K. For 
c p = 0, generalized factorization (13) reduces to usual 
representation (6). It is worth emphasizing that in view 
of the velocity independence of G on v, the modified 
operator of shifting velocities remains to be evaluated 
exactly for any b p and c p , namely, 

e^ At+c ^ At3 {r,v} = {r,v + 6 p aAt + c p GAt 3 }. (14) 

For quantum systems, where C = ^ |Vi^| 2 , the corre- 
sponding calculations also present no difficulties (at least 
for particles in external fields), because this requires only 
knowing the gradient of the potential. 

An important feature of decomposition integration 
(13) is that it, being applied to classical dynamics sim- 
ulations, conserves the symplectic map of particle's flow 
in phase space. This is so because separate shifts of po- 
sitions (7) and velocities (14) do not change the phase 
volume. The time reversibility S(-t)JZ(t) = TZ(0) of 
solutions (following from the property 5" 1 (t) = S(— t) 
of evolution operator S(t) = e ct ) can be reproduced ex- 
actly as well by imposing additional constraints on the 
coefficients a p , b p , and c p . In particular, for velocity- like 
decompositions such constraints read: a\ = 0, a p+1 = 
ap-p+i, bp = bp-p+i, and c p = cp_ p +i. Then single- 
exponential subpropagators will enter symmetrically into 
the decompositions, providing automatically the required 
reversibility. The case when the operators of shifting ve- 
locity and position are replaced by each other in the re- 
sulting symmetrical decomposition is also possible. This 
leads to a position-like integration which can be repro- 



duced from Eq. (13) at a p = ap_ p +i, b p = ap_ p , and 
c p = cp-p at bp = and cp = 0. 

The above symmetry will result in its turn to auto- 
matic disappearing all even-order terms in the error func- 
tion 0(At K+1 ). For this reason, the order K of time- 
reversible (self-adjoint) algorithms may accept only even 
numbers (K = 2,4,6,...). The cancellation of the re- 
maining odd-order terms up to a given order will be pro- 
vided by fulfilling a set of basic conditions for a p , b p , and 

c p . For example, the condition J2 P =i a p = J2 P =i bp = 1 
is required to cancel the first-order truncation uncertain- 
ties. Then the error function can be cast in the form 



0(At K+1 ) = 3 At 3 + 5 At 5 + 7 At 7 + 
... + K+1 At K+1 . 



(15) 



In order to kill higher odd-order truncation terms in 
Eq. (15), let us write down explicit expressions for 3 , 
5 , and 07 (this will be enough to derive algorithms 
up to the eighth order). Expanding both the sides of 
Eq. (13) into Taylor's series, and collecting the terms 
with the same powers of At one finds: 



O a = a[A, [A,B]] + (3[B, [A,B]} 



(16) 
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5 = 7i [A [A, [A, [AM\\+12[A, [A, [B, [A, B}}}} + l3 [B , [A, [A, [A,B]}}} + 74 [ZJ, [B, [A, [A,B]]]] , (17) 



7 = Ci[B, [B, [A, [B, [A, [B,A]}]}]} + ( 2 [B, [B, [B, [A, [A, [B,A]}]}]} + ( 3 [B, [B, [A, [A, [A, [B,A]}]}]} + 

Q[B, [A, [B, [A, [A, [B,A\]}]}] + ( 5 [A, [B, [B, [A, [A, [B,A\]}]}] + ( 6 [A, [B, [A, [B, [A, [B,A]}]}]} + (18) 
( 7 {B, [A, [A, [A, [A, [B, A}}}}}} + &[A, [B, [A, [A, [A, [B, A}}}}}} + ( 9 [A, [A, [B, [A, [A, [B, A}}}}}} + 
Cio[A,[A,[A,[A,[A,[B,A]}]}}}. 



Here we take into account the fact that the operators 
B and C commute between themselves, i.e. [B,C] = 0, 
so that any occurrence of constructions containing the 
chain [B, [B, [A, B]]] has been ignored (in particular for 
fifth-order truncation term C 5 this has allowed us to ex- 
clude the two zero- valued commutators [B, [B, [B, [A, B]]]) 
and [A, [B, [B, [A,B]]}]). The multipliers a, /3, 7i_ 4 , and 
Ci-io, arising in Eq. (16)-(18), are functions of the coeffi- 
cients a p , bp, and c p , where p — 1, 2, . . . , P. The concrete 
form of these functions will depend on P and the version 
(velocity or position) under consideration. 

The most simple way to obtain explicit expressions for 
the multipliers consists in the following. First, since we 
are dealing with self-adjoint schemes, the total number 
of single-exponential operators (stages) in Eq. (13) is ac- 
tually equal to S = 2P— 1, i.e. it accepts only odd values 
(mention that one of the boundary set of coefficients is 
set to zero, a\ — or bp = cp = 0). Then we can al- 
ways choose a central single-exponential operator, and 
further consecutively applying P — 1 times the two types 
of symmetric transformation 

e W <,, + 1) +0 < ™+ 1) _ e Aa,( n) At e W (n) +0^ At 

(19) 

e>v <™+i) +C i<"+i) _ 

c Bfc<") At+Cc (n > At 3 C W<' ,) + e> ( ™> gSb*™' At+Cc (n ~>At 3 

come to factorization (13), where 

W = (vA + aB)At 

and O is defined by Eq. (15). The quantities aS n \ 
and c(") are related to a p , b p , and c p , respectively 
(the relationship between n and p is determined below). 



For velocity-like decomposition with even P or position- 
like at odd P, the central operator is correspondingly 

c Aa (P _ 2)/2+1 At or e Aa (P _ 1)/2+1 At^ g Q tn&t nere WQ must 

put <t(°) = as well as = = j{% = = 

and let either = a( P _ 2 )/2+i or y = a (p-i)/2+i 
on the very beginning (n = 0) of the recursive proce- 
dure. The start of the procedure should be performed 
with the second line of Eq. (19) at = 6(p_2)/2 and 
c( 0) = C(p_ 2 )/2 or M 0) = 6(p_i)/2 and c^ 0) = C(i_ 2 )/2 with 
further decreasing the index p with increasing the number 
n = 1, 2, ... P - 1 at o(") = a p , b^ = b p , and = c p 
in both the lines of transformation (19). For velocity- like 
decomposition with odd P or position-like at even P, 
the central operator will be e Bb (P -i )/2+ iAt+Cc (P _ 1)/2+1 At 3 
or c Bb (P _ 2)/2+1 At+Cc (P _ 2)/2+1 At 3 ^ correS p 0n di n g to cr (0) = 
b (p-i)/2+i and (3^ = c (P _ 1)/2+1 or = b {P _ 2)/2+1 
and (3^ = C( P _2)/2+i, respectively, with = and 

= = Ci-io = 0- In this case, the procedure 

should be started with the first type of transformation at 
a(°) = 6(p_i)/2+i or = 6( P _ 2 )/2+i with decreasing p 
at increasing n for = b p , = c p , and a(") = a p in 
Eq. (19). 

The recursive relations between the multipliers u, a, 
a, (3, and 71-4 corresponding to the first line of Eq. (19) 
are: 

v (n+l) = v (n) + 2ffl („) i ^(n+1) = ^(n) ) ^0) 

a ("+l) = a (") _ flW^W ( a (») + yW) /6 , (21) 
£(»+!) =/3 (n) _ a (™) CT (™) 2 /6, (22) 



7 (n+l) = 7 (n) + o („) + ^(n)) (( 7(fl (n) 2 + 7a («) Jy (") + ^(™) 2 ) (T («) _ 6 0a< n >)/360 , 

^(n+1) = 7 (n) + a (n)( 3()a (n) (T (n) _ 30a («)^(n) _ 3 Q^)„(n) + 3 (a (n) ^(n) 2 + 2o (n) ^(n) ^(n) 2 + ^(n) 2 (T (n) 2 ^ /lgQ ^ 

(23) 

7s" +1) - 73 n) + « (n) ^ (n) ((8(a (n)2 + 12a< n >i/< n > + i^(") 2 ) ( T (n) - 120a("))/360 , 

7 (n+l) = 7 (n) + a (n )(J {n) ^(n) + Iy (n)^( n ) 2 _ g^n)) /lg0 _ 



•5 



For the second type of transformation the relations read: 



o>+i) =a (n) + fe(™) jy («) 2 / 6j 



(25) 



v (n+1) =v {n) , a {n+1) = cr (n) + 2b (n) , (24) (3 {n+1) = (3 {n) + (l2c (n) + b {n) v {n) (b (n) + cr (n) ))/6 , (26) 

7 (n+l) =7 (")_ 6 (n) iy (n) 4 / 360 , 



7 (n+l) = ^(n) _ ^(n)( 60a (n) 6 (n) _ ^(«) ( 30c (») - ftWj/W (6&(»> + <rM)))/180 , 

7s" +1) = 73° +& (n) ^")(60a^ +^") 2 (46(") - ff <"'))/360, 



(27) 



7 |" +1 ) = 7 f ) - (30aW&( n )(6W + v™) - i/W (30/3< n W n > + m^c™ - 

3&(") V n ) + 30c (rl V (n) - 26 (n) V n) er (n) - b {n) v (n) <r (n)2 ))/l%Q . 



The relations for Ci-io are presented in Appendix. In 
such a way, at the end of the recursive process (i.e. after 
P — 1 steps) the multipliers can readily be obtained. The 
form of the first two multipliers are particularly simple 



and look as v 



and a 



Ep=i V 



So that, as 



was already mentioned above, putting v = 1 and a = 1 
will cancel the first-order truncation uncertainties (be- 
cause the resulting exponential propagator must behave 
like e (- /t+8 ) A *). Next multipliers should be set to zero and 
we come to the necessity of solving a system of non-linear 
equations (so-called order conditions) with respect to a p , 
bp, and c p . We shall now consider actual self-adjoint al- 
gorithms of orders K = 2, 4, 6, and 8. 



3. Force- gradient algorithms of order two 

Putting P = 2 in Eq. (13) with oi = 0, h = b 2 = 1/2, 
a 2 = 1, and c\ = c 2 = £ leads to the following velocity- 
force-gradient algorithm of the second (K = 2) order, 



^A+B)At+0 3 At 3 



(28) 

with a = 1/12, and = 1/24 + 2£. Note that here and 
below, for reducing the number of unknowns, we will al- 
ways take into account in advance the symmetry of co- 
efficients a p , bp, and c p as well as the fulfilling the first- 
order conditions v — ^ p=1 a p = 1 = J2 P =i b p = 1 = a 
when writing decomposition formulas. Then solving the 
equation (3 = yields £ = —1/48 and we come to the al- 
ready found integrator (12). It is worth remarking that 
negative values of quantities c p at force gradients have 
nothing to do with the above problem of positiveness 
of time coefficients arising at velocities and forces, i.e., 
for a p and b p . The reason is that the incremental ve- 
locity b p &At + c p GAt 3 in Eq. (14) can be rewritten as 



(b p SL + c p GAt A 



}At = b n &At, and thus treated as the ve- 



locity changing in a modified step-size-dependent accel- 
eration field a = a + ^GAt 2 . 



The position counterpart of Eq. (28) is obtained from 
Eq. (13) at ai = 0, oi = a 2 = 1/2, b x = 1, b 2 = 0, a = £ 
and c 2 = 0, that yields 



c (A+B)At+0 3 At 3 _ e A4± e BAt+£CAt 3 e A^ 



(29) 



for which a = -1/24 and (3 = -1/12 + £. Letting 
£ = 1/12 will minimize the third-order truncation errors 
to the value a[A, [A, B]]At 3 which is even twice smaller 
in magnitude than that of the velocity version. Note, 
however, that for both versions (28) and (29), which re- 
quire one force plus one force-gradient evaluations per 
time step, the order of integration is not increased with 
respect to the usual (when £ = 0) Verlet integrators re- 
quiring only one force recalculation. In view of this, the 
applying force gradients in a particular case of P = 2 can 
be justified only for strongly interacting systems when 
the kinetic part A of the Liouville operator C is much 
smaller than the potential part B, i.e., when C = eA + B 
with e « 1. Then the remaining part a[A, [A, B}] At 3 
of local uncertainties will behave like oc e 2 and can be 
neglected. 



4- Force- gradient algorithms of order four 



Further increasing P on unity allows us to kill exactly 
both the multipliers a and (3, that is needed for obtain- 
ing fourth-order (K = 4) integrators. So that choosing 
P = 3 leads to the velocity-like propagation 



(A+B)At+0(At !i ) 



(30) 



e 2 C " 



e 2 e 



following from Eq. (13) at a\ — 0, a 2 = a 3 = 1/2, 
bi = b 3 = A, b 2 = 1 - 2A, ci = c 3 = £, and c 2 = X- 
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Here relations (21), (22), (25), and (26) come to the two 
order conditions 

1-6A 1 A A 2 

* = = 0, P = - V2+ - 2 - Y + 2Z + X = 0, 

with three unknowns A, £, and x- The first unknown is 
immediately obtained satisfying the first condition, 



.,{A+B)At+0(Af) 



(36) 



and is obtained from Eq. (13) at P = 3 with a\ = 03 = A, 
02 = (1 — 2A), bi — b 2 — 1/2, ci = c 2 = £, an d 
63 = C3 = 0. Here, the number of unknowns coincides 
with the number of the order conditions 



A 



(31) 



1 A A 2 
12-2 + T = °- 



The second equality is then reduces to 2£ + x — 1/72, 
resulting in a whole family of velocity-force-gradient al- 
gorithms of the fourth order. In general, such algorithms 
will require two force and two force-gradient recalcula- 
tions per time step. 

Remembering that we are interested in the derivation 
of most efficient integrators, three cases deserve to be 
considered. The two of them are aimed to reduce the 
number of force-gradient recalculations from two to one. 
This is possible by choosing either 



£ = 0, X = 



1 

72 



X = °' ^=144- 
In the third case we will try to minimize the norm 



7 



11 



ll 



ll 



ll 



(32) 



(33) 



(34) 



of fifth-order truncation errors 0(At 5 ) at £ 7^ and 
X = 1/72 — 2£ 7^ 0, treating £ as a free parameter. In view 
of recursive relations (23) and (27) , explicit expressions 
for the components of 0(At 5 ) are 



7 - 30A 1 
7i= r „„» , 72 : 



X A 2 £ 



1 A A 2 



5760 



480 24 24 6' 73 360 48 24' 



74 



1 
120 



A_ 

16 



™!_*! + i_x/l_A 

48 8 6 2 V 3 



Then taking into account Eq. (31) one finds the function 



solving of which yields two solutions, 



A 



it 



1 



(37) 



Then the norm of truncation uncertainties 0(At 5 ) ap- 
pearing in Eq. (36) is 7 = (1873 T 40\/2T87) 1 / 2 /2160, 
so that the preference should be given to sign "-" in 
Eq. (37), because this leads to a smaller value, 7_ w 
0.000715, of 7 (whereas 7+ » 0.0283). Position algo- 
rithm (36) needs, like velocity version (35), in two force 
and the same number of force-gradient evaluations per 
time step. 

Integrators (32) and (37) have been previously derived 
by Suzuki [B4j b ased on McLachlan's method of small 
perturbation |p3| and referred by Chin |2(| to schemes A 
and B, respectively. Algorithms (33) and (35) are new 
and will labeled by us as A' and A". While scheme A' 
seems has no advantages over the A-integrator, the new 
algorithm A" corresponds to the best accuracy of the in- 
tegration, because it minimizes 7. It requires, however, 
one extra force-gradient evaluation and, thus, can be rec- 
ommended for situations when this evaluation does not 
present significant difficulties. 

With the aim of considerable decreasing the truncation 
errors in a little additional computational efforts, Chin 
p6f has proposed to consider extended force-gradient al- 
gorithms of the fourth order. This has been achieved 
by increasing the number of force recalculations on unity 
with respect to the necessary minimum, i.e. choosing 
rtf = 3. At the same time, the number of force-gradient 
evaluations was fixed to its minimal value n g = 1. Within 
our general approach, it is possible to introduce two 
fourth-order schemes satisfying the above requirements. 
The schemes are 



7 



135 v / 2048 
with the minimum 7 m 



VT9 + 12240£ + 6480000£ 2 



17 



18000 



/66T/43200 : 
71 



0.000595 at 



X 



4500 



(35) 



At the same time, the values of 7 corresponding to first 
two algorithms (32) and (33) constitute >/l9/2048/135 « 
0.000713 and 7^/8640 « 0.00334, respectively. 
Position version of (30) reads 



c (^+B)At+C(Ai 5 ) _ A9At BXAt A(l-29)^ 



XC 



B(l-2\)At+ x CAt 3 p A(l-29)^ e B\At p A9At 



and 



e (A+B)At+0(At 5 ) _ e B\At+£CAt 3 e A9At e B(l-2\)A± 
x e A{l-29)At e B(l-2\) At & A9At e B\At+£CAt 3 



(38) 



(39) 



following from Eq. (13) at P — 4 and corresponding to 
position- and velocity-like integration, respectively. Note 
that further we will not present the relationship between 



7 



the coefficients a p , b p , c p of Eq. (13) and reduced vari- 
ables (such as, for example, 6, A, x m Eq. (38)) in view 
of its evidence. 

The order conditions for scheme (38) are 



24 I 4 



= 0. 



y 12 2 2 V J 
and solving them one obtains 



X = 



1 ± %/6A(l - A) 



12 



(40) 



Relations (40) constitute a family of extended force- 
gradient position algorithms (38) of the fourth order with 
A being a free parameter. Chin p6| has introduced an 
algorithm like (38) in somewhat another way, namely, as 
a symmetric product of two third-order schemes. This 
results only in one set of time coefficients which can be 
reproduced (at sign "-") from Eq. (40) as a particular 
case corresponding to 



X 



— (41) 
192 



and has been referred to scheme C. 

Solution (41) may not, however, be necessarily optimal 
in view of the fact that it does not minimize the norm 
7 (see Eq. (34)) of truncation uncertainties 0(At 5 ). In- 
deed, the components of 7 for scheme (38) are 



1 



7i 



73 



74 



1920 
1 



6912A 
5 



72 



'3 

:!(»0V2 vJXiA 




6±5V6A 
2880 



30A' 



where Eq. (40) has been used to express the function 7(A) 
in terms of one parameter A exclusively. The global min- 
imum of this function is 7 m i n ~ 0.000141 and achieved 
(at sign minus) at 

A = 0.2470939580390842E+00 , and thus 

9 = 0.8935804763220157E-01 , (42) 

X = 0.6938106540706989E-02 

(all results found numerically will be presented within 
sixteen significant digits for schemes up to the eighth or- 
der and within thirty two digits for order ten and higher) . 
On the other hand, the value of 7 corresponding to 
scheme C (Eq. (41)), is equal only to V87817/414720 « 
0.000715, i.e. it is approximately in 5 times larger than 
that of the optimized algorithm (42). The last algorithm 
we will designate as scheme C. 



A similar pattern is observed in the case of extended 
velocity-force-gradient integration (39). Previously Chin 
and Chen |3l"| ] have indicated that for quantum mechanics 
simulations the integration of such a type is more prefer- 
able than position-like scheme (38), because it requires a 
fewer number of spatial Fourier transforms. Again using 
the symmetric product of two third-order integrators to 
increase the order from three to four, they have obtained 
the following set 



A = 



1 



1 

384 



(43) 



of time coefficients and referred it to scheme D. We have 
realized that this set is not only possible and found a 
whole family of solutions (which includes (43)), namely, 



A= — ((;■ 
12 



1 



^=-2^ 6 



1 



1) 



where 9 should be considered as a free parameter. The 
optimal solution, which minimizes the norm 7 of fifth- 
order errors to the value 7 min w 0.000855, is 



A = 0.4432204907934768E-01 , 
9 = 0.2409202729169543E+00 , 
£ = 0.4179297897540420E-02 



(44) 



and will be labeled as scheme D'. At the same time, the 
norm of errors corresponding to scheme D (Eq. (43)) is 
equal to 7 = V237457/414720 w 0.00117, i.e. it exceeds 
the minimum, that may results in decreasing the preci- 
sion of the calculations. 

As can be ensured readily, the time coefficients arising 
at basic operators A and B under exponentials are posi- 
tive for all the fourth-order force-gradient algorithms de- 
scribed in this subsection. Therefore, contrary to usual 
force Forest-Ruth-like schemes, such algorithms can sim- 
ulate dynamical processes in all areas of physics and 
chemistry without any principal restrictions. 



5. Force- gradient algorithms of order six 

Beginning from P = 5, the force-gradient factorization 
being written in velocity representation allows to elim- 
inate the components of truncation uncertainties up to 
the sixth order (K = 6) inclusively. In view of Eq. (13), 
such a representation reads 



^(A+B)At+0(At 7 ) 



e mAt+ l j,CAf' e A8At e B\At+tCAt i 



^A(1-20)A1 8(l-2(\+-ff)) At+xCAt 3 A{l-26 



(45) 



BXAt+^CAt A8At BSAt+fiCAt 

x e e e 



The number of unknowns in propagation (45) is the same 
as the number of order conditions which now take the 
form 



8 



(3 = x - ^ ( 1 - - 6A 2 (20 - 1) - 61? + 6i? 2 - 6A(20 - l)(2i? - 1) - 24£ ) = , 



7i = 



5760 



7 - 3OA(20 - 1) 2 (1 + 46* - AO 2 ) - 30i? = 



" " = 480 1 1 ~ 2 ° X + 80fi ~ ^ ~ 86> + 18 ~ 12 -* ~ + 
2OA(20 - 1)(0 + 2i? - 60i?) + 80£ - 48O0£ + 48O0 2 £ ) = , 



(46) 



73 = t^- ( 2 - 3OA 2 (20 - l) 3 - 15i? + 30i? 2 - 15A(26> - 1) 2 (1 - 4i? - 0(4i? - 2)) ) = , 



74 = — \2 + 40m - 3OA J (20 - l) 2 - 40%(1 + A(60 - 3) - 3i?) - 15i? + 35i? 2 - 30i? J - 5A 2 (20 - 1) 

x (7 - 18i? + 60(2i? - 1)) + 5A(20 - 1)(3 - 14i? + 18i? 2 + 20(1 - 6i? + 6i? 2 )) + 40£ - 24O0£ + 48001?^ = . 



The unique real solution to system (46) is 

Q _ 1 y/675 + 75V6 5 

2 30 2-^675 + 75^' 



A = - 



50 



50 2 



£ = - — — I (471 

144 36 288 ' v ' 



X : 



1 _0_/0 
144 ~ 36 V 2 



1 



fj, = 0. 



Solution (47) constitutes a velocity-force-gradient algo- 
rithm of the sixth order with four force and three (since 
/i = 0) force-gradient evaluations per time step, i.e., with 
m = 4 and m g = 3. Its advantage over usual sixth-order 
integrators consists in the fact that it is composed from 
a considerably smaller number, namely S = 2P — 1 = 9, 
instead of 15, of single exponential operators. The norm 



10 



Ed 

k=i 



(48) 



of seventh-order truncation errors 0(At 7 ) (see Eq. (45)), 
corresponding to solution (47), is equal to ( w 0.00150. 
Note also that the position version of decomposition (45) 
does not exist at P = 5, because then the number of 
unknowns is less than the number of order equations, 
resulting in the absence of solutions. 



As has been shown in the preceding subsection for the 
case of fourth-order integration, algorithms with minimal 
numbers m of force evaluation may not lead to optimal 
solutions. The reason is that slight increasing m may sig- 
nificantly decrease the local errors and thus overcompen- 
sate an increased computational efforts. So that increas- 
ing m as well as P on unity (note that m = P — 1) and 
do not changing n g , i.e. choosing P = 6 with m = 5 and 
n g = 3, it is possible to derive from decomposition (13) 
up four (two velocity- and two position-like) extended 
sixth-order schemes. They are 



CACABABACAC . 
ABAC AC AC ABA . 



CABACACABAC , 
ACABAC ABACA , 



(49) 



where we have used an abbreviation that A and B de- 

3 BbiA *, respec- 

^BbiAt+CaAt 3 



note exponential operators e Aai and 
tively, whereas letter C corresponds to e^ 
Each of these extended schemes has itself correspond- 
ingly six, eight, four, and two sets of real solutions for 
time coefficients. We have realized that the smallest val- 
ues of the norm ( (see Eq. (48)) of local errors 0(At 7 ) 
within the sets are 0.0000264, 0.0000147, 0.000146, and 
0.00000607, respectively. So that the last scheme should 
be considered as the best. More explicit form for it is 



^A+B)At+0(At 7 ) 



e ApAt e BdAt+nCAt i c A8At e BXAt 



xc A(l-2(8+p))¥ c B(l-2(\+0))At+ x CAt 3 ^ 
xc A(i-2(8+ p )) £f c B\At e AeAt e B$At+nCAt 3 e ApAt 
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with the optimal solution 

p = 0.1097059723948682E+00 

= 0.4140632267310831E+00 

d = 0.2693315848935301E+00 , N 

51 

A = 0.1131980348651556E+01 
X =-0.1324638643416052E-01 
H = 0.8642161339706166E-03 

corresponding to £ = 0.00000607. In such a way, the 
error function has been reduced more than in 200 times 
with respect to scheme (47) for which ( w 0.00150. 

6. Force- gradient algorithms of order eight 

In the case when K = 8 we must satisfy up eigh- 
teen order conditions, namely, v = 1, a = 1, a = 0, 
(3 = 0, 7i_4 = 0, and Ci-io = 0. Taking into ac- 
count the symmetry of time coefficients a p , b p , and c p , 
this can be achieved at least at P = 12, i.e., using 
S = 2P — 1 = 23 single exponential operators. For 
P = 12 the velocity- and position- like force-gradient de- 
composition (13) transforms into the schemes 

CACACACACACACACACACACAC (52) 

and 

ACACACACACACACACACACACA , (53) 

respectively. The number of unknowns for both the 
schemes are also eighteen and we can try to solve the sys- 
tem of order conditions with respect to these unknowns. 

It is worth remarking such a system appears to be very 
cumbersome for schemes under consideration. For in- 
stance, the resulting non-linear equations of this system 
being written explicitly in Mathcmatica create a file of 
0.5 Mb in length! In view of this, our attempts to solve 
the equations symbolically have not meet with much suc- 
cess. We mention that all the results presented above for 
algorithms of orders 2, 4, and 6 have been solved analyt- 
ically or in quadratures. Saying in quadratures we mean 
that the problem was reduced to finding real zeros for 
a one-dimensional polynomial of a given order. So that 
we could identify exactly the number of solutions and 
their locations. Here the situation is somewhat different 
because we must solve the system using purely numeri- 
cal approaches, such as the Newton method. As a result, 
one cannot guarantee that we will found all possible solu- 
tions. However, solving the system on a computer during 
significantly long time, one can stay with a great proba- 
bility that we have found almost all physically interesting 
solutions and chosen among them nearly optimal sets. 

The numerical calculations has been performed in For- 
tran using the well-recognized Newton solver with numer- 
ical determination of partial derivatives. The values for 
non-linear functions (that constitute the system of equa- 
tions) were obtained using recursive relations (20)-(27), 



(Al), and (A2), but not explicit expressions for them to 
save the processor time and increase the precision of the 
computations. The initial guess for solutions were gener- 
ated at random within the interval [—2.5, 2.5] in each the 
eighteenth directions. If Newton's iterations become to 
diverge at a particular guess or during the calculations, a 
next random point was involved to repeat the process. In 
such a way, after several days of continuous attacking the 
systems of equations on an Origin 3800 workstation, wc 
found two and five solutions for schemes (52) and (53), 
respectively. The optimal among them are following 












bi -■ 


= 




0. 1839699354244402E + 00 


Ci : 


— Cl2 













6922517172738832E + 00 


b 2 = 


= 611 




0. 708438975 7230299E + 00 


C2 = 


= Cn 




0.3976209968238716E - 01 


a-3 = 


= an 




-0.3183450347119991E+00 


b 3 = 


= bio 




. 1 98 1440445033534E +00 


C i ■■ 


= ClO 




0.2245403440322733E - 01 


(24 = 


= aio 




0.6766724088765565E+00 


b 4 


= b 9 




-0.6409380745116974E-01 


Ci 


= Cg 




0.9405266232181224E - 03 




= ag 




- . 7207972470858706E + 00 


h 


= b 8 




-0.6887429532761409E + 00 


C5 


= c 8 




-0.7336500519635302E - 01 




= as 




0.3580316862350045E+00 


be 


= 07 




0.1622838050764871E+00 


C6 


= c 7 




0.2225664796363730E-01 




a 7 




-0.3756270611751488E+00 



for velocity- like integration (52), and 





bi2 









Cl2 







ai = 


= ai2 




0.41009674738801111928784693005080E+00 


bl : 


= b n 




0.48249309817414952912695842664785E-02 


Cl = 


= Cn 




0. 14743936907797528364717244 760736E - 03 


a 2 = 


= an 




-0.34123345756052780489101697378499E+00 


02 = 


= bio 




0.17492394861090375603419001374207E+00 


C2 ; 


= ClO 




0.23288450531932545357194967600155E - 03 


a 3 = 


= aio 




0.25644714021068150492361761631743E+00 


h 


= b 9 




0.29304366370957066164364546204288E+00 


C3 


= Cg 




0.61648659635535962497705619884752E - 02 


ai 


= ag 




0.27765273975812438394100476242641E+00 


h 


= b 8 




0.47448940168459770284238136482511E-01 


Ci 


= cs 




-0. 12307516860831240716732016960034E - 01 


as 


= as 




-0.56926266869753773902939657321159E+00 


bs 


= b 7 




-0.15299863411743974499219652320477E-02 


C5 


= ^ 




-0.73296648559126385387017161643798E-04 


a6 


= 0,7 




0.46629949890124853576794423820194E+00 




be 




-0.37422994259002571606842462603791E - 01 




C6 




0. 15295860994523744731993293847001E - 01 



for its position- like counterpart (53). The number of force 
evaluations per times step for schemes (52) and (53) is 
m = P — 1 = 11, whereas the number of force-gradient 
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recalculations consists n g = 10 (since c% = and thus 
the two boundary letters C in formula (52) should be 
actually replaced by B) and n g = 11, respectively. 

In view of a complicated structure of the ninth-order 
truncation uncertainties 0(At 9 ), the optimal solutions 
just presented have been chosen in somewhat other way 
than above, namely, by providing a minimum for the 
function S = ma,Xp =1 (\a p \, \ b p \). This simplified criterion 
was used, in particular, by Kahan and Li P|22j], when 
optimizing usual force algorithms. As a result, we have 
obtained S m i n = \as\ = |og| « 0.721 for scheme (52) 
and i5 m i n = |dg j = |<zg| w 0.569 for scheme (53). Since 
(5 m i n is smaller in the last case, the position-like integra- 
tion should be considered as more preferable. Its time 
coefficients have been presented even with thirty second 
significant digits to be used in applications for very ac- 
curate integration. In order to ensure that all the digits 
shown are correct in both the cases, we have carried out a 
few additional Newton's iterations in Maple with up 200 
digits during the internal computations, and taking as 
initial guesses the solutions already obtained in Fortran. 

The position-like decomposition (53) has also another 
advantage over the velocity version (52) in that the all the 
intermediate (q < P) states in velocity and position space 
stay during the integration within a given interval [0, At], 
i-e-, < J2Ui a P < 1 and < £* =1 b p < 1. This prop- 
erty may be important when solving ordinary differential 
equations (for specific physical systems or in pure math- 
ematics applications) with singularities beyond the inter- 
val of the integration (note that, in particular, any system 
of differential equations of the form d 2 x/dt 2 = f (x) is re- 



duced to the equations of motion under consideration in 
this paper). Note also that in order to construct eighth- 
order schemes within usual decomposition approach (6) 
(i.e., without involving any force-gradients), it could be 
necessary to apply up 2 • 18 — 1 = 35 (instead of 23) 
single exponential propagators. Such schemes has never 
been derived by decomposition (6) because of the serious 
technical difficulties. They can be explicitly introduced 
only by compositions of lower-order integrators (see the 
next section). Instead, using generalized scheme (13) has 
allowed us to derive eighth-order algorithms by direct de- 
compositions for the first time (the force-gradient algo- 
rithms presented in subsection II. B. 5 for order six are 
completely new as well). 

All the decomposition algorithms obtained by us in 
subsections II. B. 3, 4-> 5, and 6 are collected below in 
Table 1. Here, the designations Err3, Err5, and Err7 
relate to the norms \J o? + /3 2 , 7, and £ of correspond- 
ingly third-, fifth-, and seventh-order truncation errors 
(see Eqs. (6), (15)-(18), (34), and (48)), whereas m and 
m denote the numbers of force and force-gradient evalu- 
ations per time step. The optimal algorithms for orders 
2, 4, 6, and 8 are labeled by G2, C, G6, and G8, respec- 
tively. Among other schemes presented for each given or- 
der, such algorithms reduce the truncation uncertainties 
to a minimum. Taking into account that this reduction is 
achieved at the same or nearly the same computational 
efforts, the optimal algorithms should be considered as 
the best not only with respect to their precision but in 
view of the overall efficiency as well (see also comments 
on this in section III). 



TABLE I. The basic decomposition force-gradient algorithms 



Algorithm 



CAC 
ACA 



Order 



Err3 



Err5 



Err7 



Equations 



Remarks Label 



1 8.3310" 
1 4.17-10" 



1.34-10" 
6.48-10" 



2.24- 10 

7.25- 10 



(28) 
(29) 



New 
New + 



G2' 
G2 



BACAB 

CABAC 

CACAC 

ACACA 

ABACABA 

CABABAC 



7.13- 
3.34- 
5.95- 
7.15- 
1.41 
8.55- 



-4 
-3 
-4 
-4 

_4(a) 
-4(b) 



6.30 
2.72 
4.83 
5.59 
1.04 
2.24 



-5 
-4 
-5 
-5 

-5(b) 



(30,32) 
(30,33) 
(30,35) 
(36,37) 
(38,41/42) 
(39,43/44) 



Refs. |2J| 

New 
New 
Refs.Jp[ 

Ref. p| /New 4 
Ref. [31 1 /New 



A 

A' 

A" 

B 
C/C 
D/D' 



BACACACAB 6 4 3 

CACABABACAC 6 5 3 

CABACACABAC 6 5 3 

ABACACACABA 6 5 3 

ACABACABACA 6 5 3 



1.50 
2.64 
1.47 
1.46 
6.07 



(45,47) 
(49) 
(49) 
(49) 

(50,51) 



New 
New 
New 
New 
New 4 



G6' 

G6" 

G6'" 

G6"" 

G6 



BACACACACACA 
XCACACACACAB 

ACACACACACAC 
x AC AC AC AC ACA 



11 10 
11 11 







(52) 
(53) 



New 



New 4 



G8' 



G8 



The best algorithm within a given order 
'The value corresponding to scheme C' 
'The value corresponding to scheme D' 
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Finally, it is worth remarking that the problem of con- 
structing algorithms with only positive coefficients a p and 
bp for orders six and higher still remains. We mention 
that for order four, this problem has been resolved (see 
subsection II. B.^) by transferring the force-gradient com- 
ponent of truncation uncertainties into the exponential 
propagators. For orders K > 6, additional higher-order 
gradients should appear under these exponentials to pro- 
vide the required positiveness. Our analysis has shown, 
however, that such high-order exponentials (besides their 
very cumbersome forms) cannot be evaluated in quadra- 
tures and need in performing implicit calculations by it- 
eration. In view of this we can come to a conclusion that 
beyond fourth order, analytically integrable decomposi- 
tion algorithms with purely positive coefficients do not 
exist. Mathematically rigorous prove of this statement 
will be considered in our further investigation and pre- 
sented elsewhere. 



C. Integration by advanced compositions 

With increasing the order of integration to ten and 
higher, the construction of algorithms by direct decom- 
positions (13) becomes to be inefficient because of a large 
number of the order conditions and time coefficients. 
However, having the already derived force-gradient in- 
tegrators of lower orders K, we can try to compose them 
as 

S Q {At) = S K (diAt) . . . S K (d P At) . . . S K (diAt) (54) 

for obtaining an algorithm of order Q > K. Then the 
composition constants d p , where p — 1, 2, . . . , P, should 
be chosen in such a way to provide the maximal pos- 
sible value of Q at a given number P > 2. Note that 
lower-order propagations Sx(d p At) enter symmetrically 
in composition (54) and their total number 2P — 1 ac- 
cepts odd values. So that if a basic integrator Sk is self- 
adjoint, the resulting algorithm Sq will be self-adjoint as 
well. The idea of using formula like (54) is not new and 
has been applied by different authors in previous investi- 
gations J||,|||||] |§. But these investigations were fo- 
cused, in fact, on the compositions of usual second-order 
(K = 2) schemes (to our knowledge, no actual calcu- 
lations of composition constants for fourth- and higher- 
order-based integrators have been reported). Although 
using the second-order-based approach allowed to intro- 
duce algorithms to the tenth order further in- 
creasing Q has led to unresolved numerical difficulties 
when finding the coefficients of the compositions. 

Usually, these difficulties are obviated with the help of 
Creutz's and Gocksch's method pp| . We mention that 
according to this method, an algorithm of order K + 2 
can be derived by the triplet concatenation 

S K+2 (At) = S K (D K At)S K ((l - 2D K )At)S K {D K At) 

(55) 



of a self-adjoint integrator of order K, where Dk = 
1/(2 - 2 1 /(^+i)). l n particular, Chin and Kidwell §| 
starting from force-gradient algorithm (41) of order four 
and repeating procedure (55) up to order 12, have in- 
dicated a visible increasing the efficiency of the compu- 
tation with respect to second-order-based schemes. In 
this approach, however, the number of force and force- 
gradient evaluations (the most time-consuming part of 
the calculations) increases too rapidly with increasing K , 
namely as 3^ K ~ 4 ^ 2 relatively to the fourth-order integra- 
tor. 

The present study is aimed to overcome the above 
problems by an explicit consideration of four-, sixth-, and 
cighth-order-based (force-gradient) algorithms within 
general composition approach (54). This results in re- 
ducing the total number of basic propagations to a min- 
imum and providing significant speeding up the integra- 
tion. The composition algorithms are derived up to the 
sixteenth order inclusively. 

1. Fourth-order based algorithms 

In the case when K — 4, the basic self-adjoint propa- 
gation is 

where X\ = A + B. Explicit form of higher-order trunca- 
tion operators A5, X 7 , Ag, *n, and so on (which was pre- 
viously found for A5 and X 7 , see Eqs. (17) and (18)) are 
not important within the composition approach. Then 
formula (54) reduces to series (n = 0, 1, . . . , P — 2) of the 
transformation 

S Q n+1) (At) = S 4 {d in) At)S Q n \At)S 4 (d in) At) (57) 

with Sq = S±(dpAt) and d^ — dp- n -i- I n view of 
Eqs. (56) and (57), the structure of resulting propaga- 
tion can be cast at each n as 

SQ ^ t ^ e yiAt+y 5 At 5 +y-At 7 +y 9 At 9 +y 11 At 11 +o{At 13 ) ^ 
with 

yi =qiX u y 5 = (72*5 , yj = 9.3*7 + g 4 [*i , *i , * s ] , 

y 9 = q 5 X 9 +q 6 [X 1 ,X u X 7 }+q 7 [X u X 1 ,X 1 ,X u X 5 ], (59) 

yu = 98*n + <?9[*i) Xi, Xg] + gio[*i, Xi,Xi, Xi, X 7 ] + 
9ii[*i, *i ; Xi, Xi, Xi,Xi, X5] + qi2[Xs, Xi, X5]. 

Comparing (56) and (58) yields values of ^-multipliers at 

n = 0, namely, q{ 0) = d P , gf = d« +1 , qf = d*+\ 

gf = 4+ 5 , and gf = d^ + \ whereas, gf = qf ] = 

q 7 Q) = qf = g$ = g[? = = 0. Expanding both the 
sides of Eq. (57) into Taylor's series with respect to At, 
one finds that values for these multipliers at n > can 
be obtained using the following recursive relations 
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01 



[ n+1) - q[ n) + 2d™ , qt +1) = q { 2 ] + ^ 



K+l 



4 n+1) = ?3 n) + 2d(n) 



K+3 
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= + 2d^ K+5 , = g W + dW (g< n) + d^) (g< n >dW* +2 - g<">)/6 , 



3 (n+l) = g („) + d („) ^(„) + d („)^ ( g (n) g (n) _ ^(n) + 7g (») d (n) + 

7^d<») a - g^V * - 7^ 2 dW X+1 - 7g^d(") K+2 )/360, 



A n+1) = q { R n) + 2d(") 



if +7 



(60) 



^o +1) = ^ + (q[ n) + d {n) ) (q[ n) %i n) - 60«W + 7< Z (" ) < Z (" ) d(") + 

7^ n) d(") 2 - g<»>>)* +2 - 7g^ 2 d(") X+3 - 7^d(") K+4 )/360, 

qft+V = q[? + d<") (g^ + dW) (42g^ )2 g| n) - ^ n)4 ^ n) - 2520g< n) - llg< n)3 g< n) d< n > + 294q[ n] 'q™ 'd< n > 
42q[ n)2 q ( ^d^ 2 + 2Mqf ] d^ 2 62q^q^S n ^ + >>* ZlqifW + 
ll g <"> V>* +1 + 42q[ n)3 d^ K+2 + Q2qt )2 d^ K+ " + ZlqfW*) /15120 , 



= + d(") (g^dW 4 - q^) (q^ + d(™) 5 )/6 • 



9l2 



Applying the above relations P — 1 times will give the 
final values of g-multipliers and thus lead to the desired 
order conditions. For instance, the first condition is very 
simple and reads: qi = dp + 2 X^Lj 1 d p = 1. This pro- 
vides = X\ and guarantees (see Eqs. (56), (58) and 
(59)) that the order of the composition scheme will be at 
least not lower than that of the basic scheme, i.e. Q > 4 
in our case. All other multipliers q 2 , 93, 94, • • • gjv should 
be consecutively set to zero, forming higher-order con- 
ditions. The total number N of the conditions depends 
on a required order Q > 4 of the composition scheme. 
In particular, at Q — 6 we must kill the term 3^5 at 
fifth-order truncation uncertainties (see Eq. (58)). Tak- 
ing into account Eq. (59), this results in two order con- 
ditions, namely, q\ = 1 and q2 = which can be satis- 
fied at P = 2. Then one obtains a system of equations, 
qi = 2di + d 2 = 0, and q\ = 2df + d 2 = 0, with respect 



to two unknowns, d\ and d 2 . The system can be solved 
analytically, and the solution is di = 1/(2 — 2 1 / 5 ) = D4 
with d 2 = 1 — 2di that coincides (at K = 4) with the 
result of triplet construction (55) . This coinciding is not 
surprising because, as can be seen easily, both approaches 
(54) and (55) are identical in a partial case when P = 2 
and Q - K = 2. 

With further increasing Q, composition approach (54) 
will lead to a more efficient integration. Indeed, choosing 
Q = 8 requires the term 3^7 in Eq. (58) should be killed 
additionally. This is achieved by putting 173 = q 4 = in 
Eq. (59), and, therefore, by solving at P = 4 a system 
of four non- linear equations, qi = 1, q 2 = 0, 53 = 0, 
and <74 = with respect to the same number of un- 
knowns di, d 2 , d3, and d±. So that the minimal num- 
ber of fourth-order integrators in the eight-order com- 
position should be 2P — 1 = 7, whereas this number is 
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equal to 3^~ K ^ 2 — 9 when triplet concatenation (55) is 
used. Expressions for the non-linear equations can read- 
ily be reproduced by applying the corresponding set of 
recursive relations (60). We will not present such ex- 
pressions explicitly, because as has been realized, the or- 
der equations do not allow to be solved analytically at 
Q — K > 4 for any K > 4. But, these equations can 
be solved in a quite efficient way numerically using the 
Newton's method. Details of the numerical calculations 
are similar to those described in subsection II. B.6. Here 
(when P = 4, K = 4, and Q = 8) we have found five 
solutions, and it seems no other real solutions exist. The 
optimal set is 



di = 0.8461211474696757E+00 
d 2 = 0.1580128458008567E+00 
d 3 = -0.1090206660543938E+01 
di = 0.1172145334546811E+01 . 



(61) 



Solution (61) simultaneously leads the smallest values for 
the maximal composition coefficient max* =1 \d p \ sa 1.172 

and the norm (q^ + ql + q 2 -) 1 / 2 ~ 0.270 of the main ninth- 
order term y(At 9 ) of truncations uncertainties. 

When deriving tenth-order composition algorithms (at 
K = 4), i.e. when Q — 10, three additional order condi- 
tions arise, 95 = 0, qe = 0, and q 7 = 0, needed to elimi- 
nate the term y(At 9 ) (see Eqs. (58) and (59)). Then we 
come in overall to 7 non-linear equations which can be 
satisfied by appropriate choosing composition constants 
dp (p = 1, 2, . . . , P) at P = 7. In this case, we have iden- 
tified more than 150 real solutions and probably there are 
somebody others (we stopped the searching after several 
days of the computations). Among the solutions found 
the optimal set looks 



di = 0.80523995769578082326628169802782 
d 2 = -0.49193105914623101022388138864143 
d 3 = 0.35449258654398460535529269988483 
d 4 = -0.69573922271140223803036463461997 
d 5 = 0.39959538030329256359349977087819 
d (i = 0.54979568601438452794128031563760 



(62) 



and d-j = 1 — 2(d\ + d 2 + d 3 + d 4 + d 5 + d & ). This set 
minimizes at once the norm (9f+99+9io + 9ii + 9i2) 1 ^ 2 of 
the main eleventh-order term y^At 11 ) of truncation er- 
rors and the quantity maxj =1 \d p \ to the values 0.00412 
and 0.843 (= \d 7 \), respectively. Here, the number of ba- 
sic propagations (stages) is 2P — 1 = 13, i.e. more than 
in two times smaller than this number ^(Q~ K )/ 2 = 27 
within triplet concatenation (55). 

In order to introduce twelfth-order algorithms, Q = 12, 
on the basis of fourth-order compositions it is necessary 
to deal with P = 12 unknowns d p to fulfill the same num- 
ber of the order conditions qi = 1, and <?2-i2=0. Here 
we have found more than 200 real solutions and perhaps 
there are somebody else. The best among them, which 
minimizes rnaxj^j \d p \ to the value 0.611 (= |di2|), is 



(63) 



di = 0.17385016093097855436061712858303 
d 2 = 0.53377479890712207949282653990842 
d 3 = 0.12130138614668307673802291966495 
d 4 = 0.29650747033807195273440032505629 
d 5 = -0.59965999857335454018482312008233 
d 6 = 0.09043581286204437145871130429094 
d 7 = -0.43979146257635806886778748138962 
d 8 = -0.30251552922346495057010240779104 
d g = 0.59895872989247982114545906953712 
dw = 0.31236416538275576151816280776696 
du = -0.59081230769647833184090443445303 

with di2 = 1 — 2 J2 P =i dp- Thus, the minimal number of 
fourth-order stages needed to compose the twelfth-order 
algorithm is 2P - 1 = 23, instead of up 3 (g ~ K)/2 = 81 
as in the case of usual triplet construction (55). 



2. Sixth- and eighth- order 



algorithms 



When K = 6 or 8, the basic propagation reads 



or 



S 6 (r) 



S S (r) 



X 1 r+X 7 T 7 +X a r a +X 11 T 11 +X 13 T li + . 



^X l T + X 9 T (, +X ll T 11 +X l3 T 13 + X^T lb + . 



(64) 



(65) 



respectively. Here, the compositions reduce to the recur- 
sive transformation 

S% +1 \At) = S 6 ^(d {n) At)S^\At)S^(d {n) At) (66) 

with Sq^ being equal to Ss(dpAt) or Ss(dpAt) and 
n = 0, 1, . . . , P — 2. The left-hand-side of expression (66) 
can again be presented at each n as a single exponential, 



S Q (At) 



p yi r+y K +i r K+1 + y K+3 r K+s +y K +& r K+ ° +y K +7T J ' 



where now 



(67) 



y\ = qiXi, yx+i = 92<%k+i, 

yK+3 = q3^K+3 + 94 [^1, Xk+i] , 
yK+5 = Q5^K+5 + 96 [^l, Xr+3\ + 

q-j[X\, X\,X\, X\, Xk+i], 

yK+7 = qsXx+7 + q^[Xi,Xi,Xx+b] + 
Qw[Xi, X\, X\,X\, Xk+3\ + 
9n [Ai, Xi,Xi,Xi,X\,X\, Xk+i] 



Recursive relations for multipliers 9i_n, corresponding 
to transformation (66), remain the same in form as in 
the case K — 2. So that we should merely to put cither 
K = 6 or K = 8 in Eq. (60) to obtain the required set of 
order conditions. 
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In view of the equivalence of Eqs. (54) and (55) at Q = 
if + 2, the first step on increasing the order of composi- 
tion scheme to Q = 8 when if — 6 or Q — 10 when if = 8 
is trivial and yields P = 2, di = 1/(2 - 2 1 /< x + 1 )) = L> x , 
and c?2 = 1 — 2di. The next steps on increasing Q to 
the higher values if + 4, if + 6, and if + 8 at if = 6 
or 8 are similar to the steps described above for if = 2. 
Namely, they lead to the necessity of solving numerically 
the system of P non-linear equations, qi = 1, q 2 = 0, . . . , 
qp = 0, with P = 4, 7, and 11, respectively. The only 
difference from the case if = 2 is that at if = 6 or 8 
and Q = if + 8, the number of equations reduces from 12 
to P = 11, because of a somewhat simplified structure 
of the last truncation operator shown in Eq. (67) with 
respect to that appearing in Eq. (59). So that below we 
will present final results only with brief comments for 
each the above cases. The best set among the solutions 
found were identified as those that minimize the quantity 
5 = maXp =1 | dp | (almost always this led to the minimiza- 
tion of the norm for the main term of truncation errors 
as well). 

For if — 6 and Q = 10 there are five solutions with 
the best set 

di = 0.88480139304442862590773863625720E+00 

d 2 = 0.11922404430206648052593264029266E+00 (68) 

d 3 = -0.10677277516805770678518370004925E+01 

with di = 1 — 2(di + d 2 + d 3 ) and <5 m j n = |cLt| = 1.127 
(within three significant digits). At if = 8 and Q = 12 
we have found again five solutions and the optimal one 
is 

di = 0.90803696667238426284572611022928E+00 

d 2 = 0.95777180465215511634906238400062E-01 (69) 

d 3 = -0.10545412798113627599734519738778E+01 

with d 4 = 1 — 2(di + d 2 + d 3 ) and (5 min = |d 4 | = 1.101. 

For if = 6 and Q = 12 there were more than 150 
solutions with the optimal set 



di = 0.64725339206305240605385248392083 
d 2 = 0.44631941526959576960102601257986 
d 3 = -0.66447133641046221008529452937721 
di = -0.58260619571844248816548809046510 
d 5 = 0.64081619589013117205634311707157 
d 6 = 0.31805596598883340430918587031701 



(70) 



and dj = 1 — 2(di + d 2 + d% + d 4 + d 5 + d 6 ) with 
<5min =\d 3 \= 0.664. 

When if — 8 and Q = 14 we have computed more 
than 150 solutions also and identified among them the 
following optimal one 



di = 0.61158201716899487377123317047417 
d 2 = 0.46763050598682150405078600842681 
d 3 = -0.63245030403272077359889720182431 
d 4 = -0.58223379020720528275072356442667 
d 5 = 0.62109852451075548059651686410928 
d 6 = 0.29686555238409826518407483052733 



(71) 



with d 7 = 1 — 2(di + d 2 + ds + c?4 + d 5 + d 6 ) an d 
<5 min = \d 3 \ = 0.632. 

At if = 6 and Q = 14 the best set, among more than 
200 solutions realized, is 



di 

d 2 = 

d 3 
di ■ 



d 8 ■- 
d 9 ■ 
dio 



■ 0.32557163066085080712970217977681 
-0.47389771786834222637653653795835 
: 0.54376649763596364670254533524499 
-0.64055411141298491334240825973418 
: 0.28139025047030322588052971757542 
: 0.56345778618405675650229011409013 
: 0.64205004597526944181678051477448 
-0.16972825772391310721875128881451 
-0.57973031669054683392549871514985 
: 0.27398580283063379870623390979762 



(72) 



with dn = 1 — 2^2 p=1 dp and <5 m j n = \d-j\ = 0.642. 

Finally for if = 8 and Q — 16 the optimal solution, 
among again more than 200 sets calculated, is 



di 

d 2 
d 3 = 
di 

d 5 = 

da - 
d 7 - 



dio = 



0.29642254891413070953312450213071 
0.55268563185301488324882994018746 
-0.58134339535533393315605544309940 
0.23403665265420481243563202333267 
-0.51788958989817055303978658827453 
-0.43983975477992920522811970527874 
-0.201 370781 50942 169957468 1 1 1993444 
0.34412872002528894622975927197416 
0.03072591760996558798895428309765 
0.48652953960727041281280535031455 



(73) 



with dn = 1 — 2J2 p=1 dp and <5 m i n = |dn| = 0.592. 

As can be seen, the quantity 5 m ; n decreases with in- 
creasing Q at any if (4, 6, and 8) considered. Besides 
the improvement of the precision of the integration, this 
leads to an extension of the region of stability of the com- 
position algorithms. Indeed, the constants d p appear in 
the compositions (see Eq. (54)) in the form of the term 
d p At (and its combinations of different orders when eval- 
uating truncation uncertainties 0(At® +1 )). Then, tak- 
ing into account that 5 = maXp =1 |d p |, the maximal value 
for the size At of the time step, at which these uncertain- 
ties do not exceed an acceptable level of precision, can 
be estimated as Ai max <~ 1 /£,„;„. This also explains a 
well correlation of S m i n with the minimum for the norm 
of truncation errors. 

Sixth- and eighth-order based compositions may 
have advantages over algorithms basing on fourth-order 
schemes especially when constructing very precise inte- 
grators with high values of Q. For instance, in order to 
derive an integrator of order Q = 16 on the basis of triplet 
concatenation (55) of a scheme of order if = 4, it is nec- 
essary to apply 3(Q~ K ^ 2 = 729 fourth-order stages. Tak- 
ing into account that each such a stage requires rif — 3 
force and n g = 1 force-gradient evaluations (see Eqs. (38) 
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or (39)), one obtains in total the numbers rif = 2187 and 
n g = 729 corresponding to a whole time step. On the 
other hand, in view of result (73), an integrator of or- 
der Q = 16 can be composed at K = 8 and P = 11 
using 2P —1 = 21 eighth-order stages for each of which 
rif = n g = 11 (see Eq. (53)). So that the overall number 
of force and force-gradient recalculations will be equal 
only to 231 that is much smaller than the above values 
2187 and 729 obtained in the case K = 4. 



III. APPLICATIONS OF FORCE-GRADIENT 
ALGORITHMS 

A. Molecular dynamics simulations 

In molecular dynamics (MD) simulations we dealt with 
a system of N = 256 identical (m = particles in- 
teracting through a Lennard- Jones potential, $(r) = 
iu[(a/r) 12 - (cr/r) 6 ] . The particles were placed in a cu- 
bic box of volume V = L 3 and periodic boundary con- 
ditions have been used to exclude the finite-size effects. 
For the same reason, the initial potential was modified 
as ip(r) = $(r) — $(r c ) + (r — r c )$'(r c ) at r < r c with 
tp(r) = for r > r c , where r c = L/2 is the cut-off ra- 
dius. Then the potential ip and its first-order derivative 
ip' = dp/dr will be continuous functions any where in r 
including the truncation point r = r c . This avoids an 
energy drift caused by the passage of particles via the 
surface of truncation sphere as well as singularities of 
tp'(r) and <p"(r) at r — r c . The simulations were car- 
ried out in a microcanonical ensemble at a reduced den- 
sity of n* = Na 3 /V = 0.845 and a reduced tempera- 
ture of T* = k B T/u = 1.7. All runs of the length in 
/ = 10 000 time steps each were started from an identical 
well equilibrated initial configuration p(0). The precision 
of the integration was measured in terms of the relative 
total energy fluctuations £ = ((E — (E)) 2 )/\ (E) |, where 

E = lY,tLi mv i 2 + lE^j^y) and ( ) denotes the 
microcanonical averaging. 

The equations of motion were integrated using force- 
gradient algorithms (30), (36), and (38) of the fourth 
order within schemes A, A', B, C, and C (see Eqs. 
(32), (33), (37), (41), and (42), respectively). For the 
purpose of comparison the integration with the help of 
usual fourth-order algorithm by Forest and Ruth (FR) 
M (which represents, in fact, triplet concatenation (55) 
of second-order Verlet scheme (8)) has been performed 
as well. The corresponding results for the total energy 
fluctuations as functions of the length I — t/ At of the 
simulations is presented in subset (a) of Fig. 1 for a typ- 
ical reduced time step of At* = Atiujma 2 ) 1 ! 2 = 0.005. 
As can be seen, schemes A, B, and C exhibit a similar 
equivalence in energy conservation. This is in agreement 
with our theoretical predictions presented in subsection 
II. B.4, where the precision of algorithms has been esti- 
mated in terms of the norm 7 (Eq. (34)) of fifth-order 



truncation errors 0(At 5 ). In particular for schemes A, 
B, and C, it has been obtained 7 « 0.000713, 0.000715, 
and 0.000715, respectively. Further, as expected, scheme 
A' (7 w 0.00334) is worse in precision and leads to val- 
ues of £ which are approximately in 0.00334/0.00071 ss 5 
times larger. Note that in microcanonical ensembles the 
total energy is an integral of motion, E(t) = E(0), so 
that within approximate MD simulations, smaller values 
of £ correspond to a better precision of the integration. It 
is worth remarking also that another integral of motion, 
namely, total momentum P = X)i!=i m i v i' ^ s conserved 
exactly within force-gradient approach (13). The reasons 
are that all velocities are updated at once (see Eq. 14) 
during each stage of decompositions and the fact that 
J2iLi ft = as well as Yli=i gi = (as can be verified 
readily using the structure of Eq. (11)). 

The best accuracy in energy conservation can be 
achieved within optimized scheme C (see Fig. 1 (a)) for 
which 7 as 0.000141. It minimizes £ to a level of ~ 10~ 5 
that is in factor 0.00071/0.000141 w 5 lower than those 
related to schemes A, B, and C. At the same time, the 
usual FR algorithm leads to the worst result £ ~ 10~ 3 . 
We see, therefore, that applying force-gradient approach 
allows to reduce unphysical energy fluctuations up in two 
orders in magnitude. Let us show now that this over- 
compensates an increased computational efforts caused 
by additional calculations of the force gradients. The 
processor time used for carrying out these calculations 
(see Eq. (11)) was nearly in 3 times larger than that 
needed for evaluations of forces itself. Further, we should 
take into account that algorithm C requires rif = 3 force 
and n g = 1 force-gradient recalculations per time step, 
whereas rif = 3 and n g = for FR scheme. As a result, 
one obtains that the size At of one step within FR propa- 
gation must be in (3+3-l)/3 = 2 times shorter than in the 
case of algorithm C for spending the same overall pro- 
cessor time within both the cases during the integration 
over a fixed time interval. Finally, in view of the fact that 
the global error and thus the function £ are proportional 
to the fourth power of At, i.e., £ ~ At 4 , one finds that, at 
the above conditions, the level of conservation of the FR 
scheme reduces from 10 -3 to £ ~ 10~ 3 /2 4 . So that rela- 
tive efficiency of scheme C with respect to the FR inte- 
grator is actually equal to (10" 3 /2 4 )/10~ 5 = 100/16 w 6. 

In order to ensure that scheme C (Eq. (42)) is indeed 
the best among whole family (40) of C'-like integrators 
(38), we carried out additional simulations in which the 
parameter A, being constant within each simulation, var- 
ied from one run to another. The total energy fluctua- 
tions obtained in such simulations at the end of the runs 
for two (fixed within each run) undimensional time steps, 
namely, At* = 0.0025 and 0.005, are shown in subset (b) 
of Fig. 1 as functions of A. As can be observed, the de- 
pendencies £ (A, At) have the global minimum located at 
the same point A « 0.247 independently on the size At 
of the time step. This point coincides completely with 
the minimum given by Eq. (42) for the function 7(A) 
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(see Eq. (34) and relations following just after Eq. (41)) 
which is also included in the subset and plotted by a 
dashed curve (where an upper lying part of the curve 
corresponds to sign plus in Eq. (40), whereas a lower ly- 
ing part as well as the simulation data are related to sign 
minus). So that our criterion on measuring the precision 



of the integration in terms of the norm of local truncated 
uncertainties is in excellent accord. Moreover, the en- 
ergy fluctuations appear to be proportional to that norm 
7 as £ (A, At) ~ 7At 4 and the coefficient of the resulting 
proportionality almost does not depend on A and At. 



10 
10 





f FR 




A' 


L C 






A 

C 





5000 

t/M 



(a) 



10000 



10 
10 
10 
10 
10 
10 







J. _ u - - " 

2 - ' " 




■ \ ^ 

3 \ 

4 ■ \ ' V ' 




> 

ID 




6 \ \j 





(b) 



0.0 0.2 0.4 0.6 0.8 1.0 

X 



FIG. 1. (a) The total energy fluctuations £ as functions of the length of the simulations performed using force-gradient 
algorithms A, A', B, C, and C' in comparison with the result of usual FR scheme, (b) The fluctuations £ obtained within 
extended C'-like scheme as depending on free parameter A at two fixed time steps. The values of A related to schemes C and 
C' are marked by the vertical lines. The function 7(A) is plotted by the dashed curve (see the text for additional explanations). 



It is worth remarking that the results reported in this 
subsection should be considered as the first attempts of 
applying force-gradient algorithms to MD simulations. 
In previous studies, algorithms of such a kind have been 



tested for classical [26 29] and quantum 131f] mechanics 



systems composed of a few bodies only (or even one body 
moving in an external field). The present investigations 
have demonstrated that force-gradient algorithms can be 
used with equal success in statistical mechanics simula- 
tions dealing with a great number of particles, i.e., when 
N ^> 1. In the last case, the calculations of force gradi- 
ents also presents no difficulties. Indeed, during the in- 
tegration we should first evaluate usual forces fi for each 
particle i, where i = 1, 2, . . . , N. This involves a number 
of operations which is proportional to the second power 
of N. Then in view of the structure of Eq. (11) and 
taking into account the fact that particle's accelerations 
b. l = fi I mi are already known quantities, the calculations 
of gradients gi will require a number of operations which 
is proportional to the same power of N, i.e., oc N 2 (but 
not to oc iV 3 , as it may look at first sight). Further re- 
ducing the computational efforts is possible taking into 
account that function g(r^ ) (see Eq. (11)) decreases with 
increasing the interparticle separation m more rapidly 
than initial potential tp(rij). In such a case, a secondary 
cut-off radius R c < r c can be introduced when evaluating 



gi (that is equivalent to putting g(i"y) = at ry > R c ) 
in order to speed up the calculations. 



B. Celestial mechanics simulations 

One of the simplest way to test force-gradient algo- 
rithms of higher orders is to apply them to solution of 
the two-dimensional Kepler problem. In particular, this 
way has been chosen by Chin and Kidwell |2^,|2^| when 
testing fourth-order algorithms A, B, and C and higher- 
order iterated counterparts of the last scheme. As has 
been established, this force-gradient scheme is particu- 
larly outstanding and appears to be much more supe- 
rior than usual non-gradient integrators, such as fourth- 
order by Forest and Ruth j| as well as by Runge and 
Kutta sixth-order by Yoshida 0, etc. In this 

subsection it will demonstrated that further significant 
improvement in the effectiveness of the integration can 
be reached replacing standard iteration procedure (55) 
by advanced composition approach (54). Moreover, us- 
ing our new sixth- and eighth-order force-gradient algo- 
rithms as the basis for the composition has allowed us to 
perform the computations with extremely high precision 
which exceeds by several orders the accuracy observing 
within standard fourth-order based schemes. 
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We will consider a motion of a particle (planet) of mass 
mi moving in the (gravitation) field (p(r) = —c/r of the 
central body (sun) with mass to 2 ^> mi, where c > 
is the constant responsible for intensity of the interac- 
tion. For simplifying the calculations, one neglects the 
influence of all other (i = 3, 4, . . . , N) particles (planets, 
for which m t -C TO2) in the (solar) system. Then the 
motion can be described by the following system of two 
equations, 



where r = ri — T2, and for clarity of the presenta- 
tion we have used units in which the reduced mass 
rxi\mil {rn,\-\-m2) and the interaction constant c are equal 
to unity. Since the quantity E = v 2 /2 — 1/r (which is 
associated with the total energy) presents an integral of 
motion for equations (74), it should be conserved during 
the integration. However, this will so if these equations 
are solved exactly. In numerical simulations, the local 
truncation uncertainties 0{Afi +1 ) accumulate step by 
step of the integration process, leading at t ^S> At to the 
global errors O(At^), where Q denotes the order of a 
self-adjoint algorithm. So that the quantity E can be 
presented as a function of time as 

E(t) =E Q +E Q (t)At® +0(At®+ 2 ), (75) 

where Eq = E(0) and Eq is the main step-size indepen- 
dent error coefficient. 

In our simulations we solved two-dimensional Kepler 
problem (74) with the same initial conditions r(0) = 
(10,0) and WO) = (0,1/10) as those used by previ- 
ous authors ]26| , |29| to make comparative analysis more 
convenient. The resulting highly eccentric (e = 0.9) 
orbit provides a nontrivial testing ground for trajec- 
tory integration. The numerical effectiveness of each al- 
gorithm was gauged in terms of main error coefficient 
E Q = Hm A t->o [-£(*) ~ E }/At® (see Eq. (75)). It can 
actually be extracted from the fraction [E(t) — Eq}/ At® 
by choosing smaller and smaller time steps At to be en- 
titled to completely ignore next higher-order corrections 
0{AtQ+ 2 ). This typically occurs in the neighborhood of 
At ~ F/5000, where P = tt/(2|£ , | 3 ) 1 / 2 is the period of 
the elliptical orbit. Since we are dealing with algorithms 
of high orders Q and small step sizes At, all the calcu- 
lations have been carried out in Fortran using quadru- 
ple (instead of double, as for MD simulations) precision 
arithmetics for ensuring the correctness of the results. 

The normalized energy deviations Eq/Eq obtained in 
the simulations applying fourth-, sixth-, eighth-, tenth-, 
twelfth-, and fourteenth-order algorithms are plotted in 
subsets (a), (b), (c), (d), (e), and (f) of Fig. 2, respec- 
tively, as functions of time t during one period P of the 
orbit. These deviations are substantial only near mid 
period when the particle is at its closest position to the 



attractive center. Note also that within symplectic inte- 
gration, the nonconservation of energy for periodic orbits 
is periodic and its averaged (over times t ^> P) value is 
bounded and independent on t (the independence of av- 
eraged energy fluctuations at t 3> At has already been 
demonstrated in MD simulations, see Fig. 1). That is 
way we presented the results in Fig. 2 within a narrow 
region of time near t ~ -P/2, where the maximal devia- 
tions of Eq will give a main contribution to the overall 
fluctuations. 

In the case of fourth-order integration we used most 
typical algorithms A, B, C, and C (see Eqs. (32), (37), 
(41), and (42), respectively). As can be seen from subset 
(a) of Fig. 2, the pattern here is somewhat different than 
that in MD simulations (please compare with Fig. 1 (a)). 
The algorithm C is clearly better than schemes A and B, 
that confirms the conclusion of Ref. p6[ . On the other 
hand, integrator C does not exhibit an improved pre- 
cision in energy conservation with respect to scheme C. 
Nearly the same was seen when iterating these algorithms 
to higher orders with the help of triplet construction (55). 
In particular, the sixth-order C'-counterpart appeared to 
be even slightly better than the corresponding counter- 
part of scheme C (see subset (b) of Fig. 2). At the same 
time, higher-order integrators basing on schemes A and 
B were definitely worse. So that the obvious candidates 
for fourth-order based iterations (55) and compositions 
(54) are schemes C and C. 

In order to understand why scheme C does not lead to 
the expected improvement over scheme C in this particu- 
lar situation, it should be taken into account that we deal 
with a small system, actually with one body moving in 
an effective external field. Moreover, such a body moves 
periodically and, thus, covers only small part of phase 
space during its displacement. This is contrary to many- 
body statistical systems, where the phase point may visit 
considerably wider regions of phase space. In the latter 
case, during the averaging along the phase trajectories, 
different components 71-4 of fifth-order local uncertain- 
ties (see Eq. (17)) will enter with approximately the same 
weights when forming the total error vector 0(At 5 ). This 
has been tentatively assumed when writing the norm 7 of 
that vector in the form of Eq. (34) and further minimizing 
7 to obtain algorithm C. In the case of a few-body sys- 
tem, especially with periodic motion, the above weights 
may differ considerably. This complicates an analysis 
of the truncation terms and makes impossible to find an 
exact global minimum for them within any analytical ap- 
proach. Note, however, that even here, the assumption 
on uniform contribution of truncation-error components 
works relatively well. Indeed, in view of dependencies 
shown in subsets (a) and (b), we can stay that both the 
schemes C and C are comparable in precision. The same 
was observed for their higher-order counterparts. For 
this reason (and to reserve more free space for other de- 
pendencies), in further subsets (c)-(f) we will draw only 
curves corresponding to scheme C. 
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FIG. 2. The normalized energy deviation of a particle in a Keplerian orbit. The results obtained within fourth-, sixth-, 
eighth-, tenth-, twelfth-, and fourteenth-order algorithms are shown in subsets (a), (b), (c), (d), (e), and (f), respectively. The 
basic algorithms used are: fourth-order schemes A, B, C, and C', as well as sixth- and eight-order integrators (correspondingly 
marked as G6 and G8). The curves related to higher-order algorithms concatenated on the basis of schemes C by standard 
iterations are labelled by the same letter C in each the sets. The fourth-, sixth-, and eighth-order based algorithms constructed 
within advanced composition approach are marked as S, G6 and G8, respectively (see the text). 
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When considering the sixth-order integration, we real- 
ized that direct velocity-like scheme defined by Eqs. (45) 
and (47) is much worse (the maximum deviation of Eq 
were more than in two orders larger) than its extended 
position- like counterpart given by Eqs. (50) and (51). 
This is in agreement with a prediction of subsection 

11. B. 5. The result corresponding to the position-like al- 
gorithm is plotted in subset (b) of Fig. 1 by the bold 
dashed curve marked as G6. As can be seen, all three 
curves shown in this subset, namely C, C, and G6 are 
close enough to each other. But algorithm G6 uses only 
rif = 5 force evaluations per time step, instead of rif = 9 
needed for iterated C- and C'-like schemes (for all these 
three cases the number of force-gradient evaluations is 
the same and equal to n g = 3). Therefore, for order 
six, direct decomposition approach (13) leads to more 
efficient integration than concatenations of fourth-order 
schemes. 

Beginning from order eight, the above concatenations 
based on standard iterations (55) and advanced compo- 
sitions (54) will result in completely different integra- 
tors. The simulation data for these iterated and com- 
posed C-based integrators are shown in subset (c) of 
Fig. 2 by thin (marked simply as C) and bold (marked 
as S) solid curves, respectively. The curves related to 
tenth- and twelfth-order iteration and composition inte- 
grators (based on the same fourth-order scheme C) are 
plotted correspondingly in subsets (d) and (e) of Fig. 2, 
and marked by the same letters C and S. We mention 
that C-marked curves have already been presented in 
the work by Chin and Kidwell |^6|,^9| up to order 12. 
They are redisplayed by us in order to illustrate the ev- 
ident superiority of our new composition approach over 
the standard iteration method. Indeed, for the iteration 
integrators (C-marked curves) of orders Q — 8, 10, and 

12, the magnitudes of the normalized energy coefficient 
Eq/Eq after one period are 1.44, 19.24, and 424.8, re- 
spectively. On the other hand, the magnitudes related 
to the composition integrators (S-marked curves) consti- 
tute correspondingly 0.0953, 0.0577, and 1.41, i.e., they 
are approximately in 15, 330, and 300 times smaller. In 
addition, the composition integrators are faster with re- 
spect to their iteration versions in factors 9/7, 27/13, and 
81/23 for Q — 8, 10, and 12, respectively (see subsection 
III.A.i), and thus the resulting efficiencies will increase 
yet. 

What about sixth- and eight-order-based composition 
schemes at Q > 8? First of all, let us consider the 
case of eight-order integration. Here, the direct scheme 
chosen was position-like integrator (53) (it leads to bet- 
ter energy conservation with respect to its velocity-like 
counterpart (52)). The result corresponding to this in- 
tegrator is plotted in subset (c) of Fig. 2 by the dashed 
curve marked as G8. As can be seen, the fourth-order- 
based composition scheme (S-curve) is better at Q = 8 
with respect to both direct G8 and iterated G6-like ver- 
sions. With increasing the order to 10 and 12, all they 
become to be nearly equivalent in the accuracy of energy 



conservation. But fourth-order-based approach requires 
somewhat fewer number of operations. For instance, for 
order 12, one obtains that the numbers of force and 
gradients evaluations per time step are equal for it to 
m = 23 • 3 = 69 and n g = 23, respectively, whereas these 
numbers for sixth- and eighth-order-based compositions 
G6 and G8 are rif = 13 • 5 = 65, n g = 13 • 3 = 39 and 
rif = rig = 7 • 11 = 77 (where G6- integrator requires less 
operations than G8-scheme). However, beginning from 
order 14, the situation reverses. The fourth-order-based 
composition S-approach becomes to be not longer acces- 
sible (because of the absence of explicit expressions for its 
time coefficients here). On the other hand, applying the 
standard fourth-order-based iteration C-method is very 
inefficient. In particular, at Q = 14 the maximal energy 
deviations within this method is \En/Eo\ max — 9901 
with n f = 21 • 5 = 729 and n g = 21 • 3 = 243. At the same 
time, the higher-order-based composition schemes lead to 
much accurate results, namely, |i?i4/-Eo|max = 2.065 with 
m = 21 • 5 = 105 and n g — 21 • 3 = 63 for G6- as well 
as \Eu/E \ max = 0.101 with n f = n g = 13 • 11 = 143 for 
G8-based schemes (where the better precision for the last 
scheme compensates to some extent its increased values 
for quantities rif and n g ). We see, therefore, that the rela- 
tive efficiencies of G6- and G8-based schemes with respect 
to C-approach constitute about 10 4 -10 5 . Finally, in the 
case Q = 16 (not shown in Fig. 2) we have obtained the 
values \E 16 /E \ max = 2.43 • 10 5 and \E W /E \ m ^ = 48.16 
corresponding to C- and G8-based schemes, respectively. 
Taking into account the numbers of rif and n g for these 
schemes presented at the end of subsection II. C.2, one 
can conclude that the efficiency increases here also ap- 
proximately in 10 4 to 10 5 times. 



IV. CONCLUDING REMARKS 

In this work we have formulated a general theory of 
construction of force-gradient algorithms for solving the 
equations of motion in classical and quantum systems. 
This has allowed us to extend considerably the class of 
analytically integrable symplcctic schemes. The new al- 
gorithms derived include self-adjoint direct decomposi- 
tion integrators of orders two, four, six, and eight as well 
as their composition counterparts up to the sixteenth or- 
der in the time step. As has been proven theoretically 
and confirmed in actual numerical simulations, these al- 
gorithms lead to significant improvement in the efficiency 
of the integration with respect to existing force-gradient 
and non-gradient schemes. It has been demonstrated 
that force-gradient algorithms can be used with equal 
success as for describing the motion in few-body classi- 
cal and quantum mechanics systems as well as for per- 
forming statistical molecular dynamics observations over 
many-particle collections. In all the cases the calcula- 
tion of force gradients presents no difficulties and requires 
computational efforts comparable with those needed to 
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evaluate usual forces itself. The new algorithms may 
be especially useful for the prediction and study of very 
subtle dynamical effects in different areas of physics and 
chemistry including the problems of astrophysical inter- 
est, whenever the precise integration of motion during 
very long times is desirable. 

The algorithms introduced exactly reproduce such im- 
portant features of classical dynamics as time reversibil- 
ity and symplecticity. This explains their excellent 
energy conservation and stability properties. In this 
context it should be mentioned another class of (non- 
gradient) integrators recently developed [Q on the basis 
of a modified Runge-Kutta approach. Like the force- 
gradient algorithms, the Runge-Kutta- like integrators 
also allow to produce time-reversible and symplectic tra- 
jectories in phase space with, in principle, arbitrary order 
in precision. However, such integrators are implicit and 
require cumbersome systems of globally coupled (via po- 
sitions and forces of all particles) nonlinear equations be 
solved by expensive iterations at each step of the inte- 
gration process. Since in practice such equations cannot 
be solved exactly, the time reversibility and symplectic- 
ity can be violated. This may lead, in particular, to in- 



stabilities in long-term energy conservation, i.e., to the 
same problem inherent in the tradition (nonsymplectic) 
Runge-Kutta method (see the Introduction). All these 
disadvantages are absent in the present approach, where 
the phase trajectories are propagated explicitly in time 
by applying consecutive simple shifts of particles in ve- 
locity and position space with exact preservation of the 
phase volume and reversibility of the generated solutions. 

The approach presented can also be adapted to the in- 
tegration of motion in more complicated systems, such 
as systems with orientational or spin degrees of freedom, 
etc., where splitting of the Liouville operator into more 
than two parts may be necessary to obtain analytically 
solvable subpropagators. These and other related prob- 
lems will be considered in a separate investigation. 
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Appendix 

The recursive relations for the highest-order multipli- the first type of self-adjoint transformations given by the 
ers Ci-io (see Eqs. (13), (15), and (18)) corresponding to first line of Eq. 19 are: 

c (n+l) = c („) + a (n)( 630/3 (n) 2 + ^O^'V") - 63/3< n > (6a< n > + V^)^ + 

a (») 3 (2i a («) + (27(a^ 2 + 9a< n M n > + v (n) 2 )a (n) ))/ 3780 > 



c („+i) = c („) + a (n)(33 6/3 H( 6a H + v (n)yn)> _ 5om (n) 2 _ 50 40 7 | n V") - 
cr( n ) 3 (336a^ + (l20(a^ 2 + Ua^u^ - v^ 2 )a^)) / '45360 , 



c („+l) = c (n) _ a (n)( 5Q40a (n)p(n) + a (n) (gQ^t") „ 84/3 (™y™) 2 + 72(a< n >V n > * + V"^ + 

24a ( ™V™ ) (*/" ) <7 ( " )2 - 42/3 (n) ) + (a (n)2 (88^ (rl) t7 (n)2 - 672/3 (n) )))/15120 , 

.(n+1) An) n {n)( m * n {n)(e n a{n) _ ( ,{n) , ,»\_(n) 2 \ , An) A nn o n ,>) , KfUfWW _ 1fW™),» 2 



Cl r ; = CP + a^(l68a^(60f3^ - (6a<"> + v^)a^ ) + a™ (10080^ + 5040^ - 1680( n V n > + 
192(a^ V"> 2 + 5*>> V"> 2 + 6a< n >i/< n > (l3*>V") 2 - 336/3^) + 
(a ( " )2 (272^ (n) a (n)2 - 1344/3 (n) )))/15120 , 



c („+i) = c („) _ a (n)(2520 7 |"VW + 7560 7 <™V" ) - 294p^u^ a^ + 180(a (n) - *»V n > + 

84a^(l20/3 ( ™ ) + (3^ (rl) - 22a ( ™ ) )a ( ™ )2 ) + {a^ 2 (234i/ r V™ )3 - 1512/3^V (rl) ) + 



6a (n) (420 7 | n) - 308/? ( "V'V" ) + 3i/ n) V n)3 ))/45360 , 



(Al) 
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c (n+i) = c (n) +a («)(i 8(a (n) 3 CT (™) 3 - 84a< n > (15/3^ - ( fl W + i/W) ff W J ) + (a< n > 2 (l5i/< n >a< n > 3 - 252/?< n V n >) 
6a ( " ) (210 7 |" ) - 28/?( n V n V n > + i/-)V B ' S ) + 2(630 7 | n V") - 
6304 n) v {n) - 42/?( n M n > V n > + z»V™) 3 ))/7560, 



c („+l) = c („) + fl( „) ( 2520a (r l )^ _ g4a (n) ( 8 ( a (n)^ + 12a W v W + v (^ a (n) + a (n) ^ 04 Q^) + 

(48(a(™) 4 + 120(a^V n ) + 92(a<") V") 2 + 18a (? V" )3 + !>) 4 )<7 (n) ))/ 15120 , 



c (n+l) = c („) _ a (n)( 504Qa (ny + ^Q^M _ ^(n) ^(n) + 25 20 7 < n) <7 (n ) - 420a (n) 0< n > (a ( " } + 2 V ^)<J^ + 

69(a (n) V n)2 + i/W 4 ffW a + 2( a (™> V™' (53z/ n V n > 2 - 294/3^) + (a( n > 3 (l48i/( n V n > 2 - 294/3 (rl) 
6aW (4207^ - 56^^ + 3z>> V") 2 ))/15120 , 

c (n+l) = c (n) + a (n)( 2520a (n) 2 _ 42a (n) (g (a (n) 2 + ^W^n) + ^nJ^W + H 4 ( fl W>) a _ 4 (a^ 3 (l47/3^ - 

59i/ n V< n > 2 ) + (a( n >V n >(l73i/( n V< n > 2 - 1176/3^) + 24a^ (210 7 < n) + 105 7 < n) - 28/3< T V™) 2 - 
2! »V") 2 ) + „M (5040 7 < n) + 2520 7 < n) - 84/3("V"> 2 + 5*»V") 2 )) / 15120 , 

= ([^ + a (") ( a W + i/< n >) (2520 7 < Il) - 42a^ (7( a <"> 2 + 7a< n >i/< n > + *> )2 ) + 

(31(a (n)4 + 62(a(") 3 ^(") + 42(a(") 2 ^(") 2 + lla ( "V n)3 + i/W 4 )^")) /15120 . 

For the transformation of the second type (see the second line of Eq. (19)) we have obtained: 

c (n+l) = c (n) _ ^ h (n)\{n)* + l bh {n)\(n)\{n) + A2c (n) v (n) ^(n) + 30c (n) _ .(n)^) 2 ) _ 

84a (n) (l5f3 (n) b (n) +30b (n) c (n) - 36 ( " )3 z/ (n) + lhc {n) a {n) - 2& ( ™ ) V™ ) cr (n) - b {n) v {n) a {n)2 ) - 
66<") 2 (210 7 < n) + */") 2 (l4/?(") + 63c<"> - ^ ( "V") 2 )) + &(™)(l260 7 i"V") - 
2ct ( " ) (630 7 ^™ ) +^ (n)2 (42/3<") +84c<") - !>V"> 2 ))))/7560 , 



c (n+l) = c (n) + ( l2b (n)\{n)* _ m (n)\(n)\{n) + A2c (n) v (n) (^Q/jH + 120c< n > - I/<"V n ) 2 ) - 252a<"> (20/3< n >6< n > + 

406(™) c («) _ 3& (n) V n) + 20c (,l) cr (,l) - 2b {n)2 v (n) a (n) - b {n) v {n) a {n)2 ) + 246 (n)2 (315 7 ^ } - 
j,(«) 2 ( 2 l/3(") +42c<") + v W a W 2 )) + M™' (2520 7 |" V™' + <t^ (7560 7 ^ n) - 
i/W 2 (294/3 (ri) + 168c< n > + v {n V n) 2 ) ) ) ) /45360 , 

c (n+l) = c (n) _ ( 2520a («) 2 fe (n) + ^(^(^ _ g4()a (n) I/ (n) (g c (n) _ b (n)\(n)j + 42c (")^(") 3 (7 (") - 12b^ 2 (210 7 ^ - 

i/WVW) - (2520 7 < ri V™ ) - 42/3(™^(™) 3 + 336c ( "V™ )3 + 2520 7 ^ ) ct ( " ) + j» V n ) 2 ))/15120 , 

c (n+l) = c |«) + ( 5040a («) 2 fe (n) _ 42a (n) I/ (n) ( 12Qc (n) _ j( B ) v (n) (gg^n) + ^(n))) + ^(n) ( 96fo («) V™) 3 + 84 C < n >I/< n > V n > 

+18&< n >V") V n > - (5040 7 * Il) + 2520 7 * ri) - v {n) ^ '\80 n) - 672c<") - 5z/ n V") 2 ))))/15120 , 

(A2) 

c („+i) = c („) _ ( 2520a (") 2 6 (n) _ 366 («) 3 I/ (n) 4 + 42c<"M") V") + 306(">V n >V n > + 168a™ (I5c< n > - 



6(«) ! y(")(6&(") +<t^)) - ftW (l5120 7 ^ n) i» - ^") J (252/3 (n) +504c^ W"V™^)))/45360 



3 , 
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C<" +1) = C { 6 n) - (630a (n)2 fe (rl) + 27& (n) V" )4 - 21c (n) *> )3 V"' + 96^ V"> V n > - 63o>V n) (20c (n) - 

&(»)„(") (6fc(") +cr (n) )) - & ( " ) (l260 7 i,™V" ) + i» 3 (21/?(") + 252c (n) - */™V n ) 2 )))/3780 , 

c („+l) = c (n) _ b (n) v (n)( 252Ql (n) _ ^(n) ^(n) 2 _ „(«) 4 ( 66 («) _ ) ) /L5120 , 

= C< n) + (50406(") 7 ^V(") - 42c<"M") 4 - 6fe("> V™' 5 + &( n M n ) V n >) /15120 , 

c (n+l) = C («)_ i ,(«) 3 ( 84a (r l ) & (n) _ ( 84c (n) _ jW^W^iW + 5^))) /15120 , 

C^ +1) =do ) -& (n) ^ ( " )6 /15120. 



All these relations as well as other symbolic expres- 
sions presented in the work has been carried out using 
Mathematica 4.0 and Maple 6 packages installed on the 



Silicon Graphics Origin 3800 workstation at Linz Uni- 
versity. The numerical calculations have also been per- 
formed there. 
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